%Math 316 - Spring 2017
%Numerical ODEs
%MATLAB code for EXAMPLE 3
%Written by M. Motamed @ UNM 
%====================================================================
%Part I - plot the exact solution   
%====================================================================
clear all

t0=0; %initial time
T=10; %final time
t=linspace(t0,T,1000);
y1exact=(cos(t)+sin(t)).*exp(-0.5*t);
y2exact=(cos(t)-sin(t)).*exp(-0.5*t);

figure
hold on
plot(t,y1exact,'b-',t,y2exact,'r-')
%legend('y1(t)','y2(t)')
xlabel('t')
ylabel('y(t)')
set(gca,'fontsize', 20); 
box



%=====================================================================
%Part II - compute and plot the numerical solution by RK1
%=====================================================================
P=[-.5 1;-1 -.5];
f=@(t,y) P*y;   %y is a 2x1 vector y=[y1;y2]. Hence f is also a vector   

h=1e-4;        %step size
t=(t0:h:T);    %vector of descrete grid points
n=round((T-t0)/h); %number of sub-intervals

   
y=zeros(2,length(t));   %to pre-allocate a vector for the solution
y(:,1)=[1;1];           %to set the initial solution at t0=0
for i=1:n
   K=f(t(i),y(:,i));
   y(:,i+1)=y(:,i)+h*K;
end
plot(t,y(1,:),'b--',t,y(2,:),'r--')

legend('exact y1', 'exact y2', 'numerical y1','numerical y2')


%=====================================================================
%Part III - compute and plot the numerical error
%=====================================================================
y1exact=(cos(t)+sin(t)).*exp(-0.5*t);
y2exact=(cos(t)-sin(t)).*exp(-0.5*t);
Error=abs(y1exact-y(1,:))+abs(y2exact-y(2,:));
figure 
semilogy(t,Error,'m')
set(gca,'fontsize', 24); 
xlabel('t')
ylabel('Error(t)')
box
box
