% Solve u'=-v, with u(0)=1
%       v'=u,  with v(0)=0
% Exact solution is u(t)=cos(t), v(t)=sin(t), so u^2+v^2=1
% Use theta method with different values of theta

N=100; % N+1 is number of timesteps used 

% first do explicit Euler
theta=0;
[t,U]=theta_solver(N,theta);
figure(1),clf
plot(U(1,:),U(2,:),'r','linewidth',2), hold on
figure(2),clf
plot(t,U(1,:).^2+U(2,:).^2,'r','linewidth',2), hold on

% now do implicit Euler
theta=1;
[t,U]=theta_solver(N,theta);
figure(1)
plot(U(1,:),U(2,:),'g','linewidth',2)
figure(2)
plot(t,U(1,:).^2+U(2,:).^2,'g','linewidth',2)

% finally Crank Nicolson
theta=0.5;
[t,U]=theta_solver(N,theta);
figure(1)
plot(U(1,:),U(2,:),'k','linewidth',2)

% label plots etc
axis equal tight
axis([-1.2,2,-1.2,2]), grid
xlabel('u','fontsize',20), ylabel('v','fontsize',20)
legend('explicit Euler','implicit Euler','Crank Nicolson','fontsize',15)
figure(2)
plot(t,U(1,:).^2+U(2,:).^2,'k','linewidth',2)
xlabel('t','fontsize',20), ylabel('u^2+v^2','fontsize',20)
legend('explicit Euler','implicit Euler','Crank Nicolson','fontsize',15,'Location','NorthWest')
grid

% function to solve the problem approximately using the theta method
function [t,U]=theta_solver(N,theta)
t=linspace(0,2*pi,N+1); 
dt=t(2)-t(1);
U=zeros(2,N+1);
U(1,1)=1;
L=[1 theta*dt; -theta*dt 1];
R=[1 -(1-theta)*dt; (1-theta)*dt 1];
for n=1:N
    U(:,n+1)=L\(R*U(:,n));
end
end