%Math 316 - Spring 2017
%Numerical ODEs
%MATLAB code for EXAMPLE 4
%Written by M. Motamed @ UNM 
%====================================================================
clear all


f=@(t,y) [y(2);y(3);-y(3)^2-y(1)*y(2)+sin(t)];   
%Note: inputs to f are t and y 
%t is scalar and y is a 3x1 vector y=[y1;y2;y3] 
%Hence f is also a 3x1 vector   


t0=0;               %initial time
T=2;                %final time
h=0.001;             %step size
t=(t0:h:T);         %vector of descrete grid points
n=round((T-t0)/h);  %number of sub-intervals

y=zeros(3,length(t));    %to pre-allocate a vector for the solution
y(:,1)=[1;0;2];          %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

figure
plot(t,y(1,:),'b-',t,y(2,:),'r-',t,y(3,:),'k-')
legend('y1(t)','y2(t)','y3(t)')
xlabel('t')
ylabel('y(t)')
set(gca,'fontsize', 22); 