Break up the following file into two functions and a script, saving the
functions under the appropriate file names.

_________________________

function z=fa(y,t)
global a w
z(1)=y(2);
z(2)=-y(1)+a*cos(w*t);

_________________________

function y=yexact(t)
global a w
y=-a/(w^2-1)*(cos(w*t)-cos(t))+cos(t);

_________________________


clear
echo on
%This program demonstrates the use of the Euler
%method, backward Euler method, trapezoidal rule,
%2-stage runge kutta, and 4-stage runge kutta rules
%in the integration of an arbitrary function f(y,t).
%We will also Newton's method to solve the implicit 
%parts of the
%trapezoidal rule and backward Euler methods.
pause
%The specific function we will integrate is the undamped
%forced oscillator y"=-y + a*cos(w*t).  The forcing
%function is parameterized by the two parameters a
%and w.  In order to pass these parameters to 
%the matlab function, we will simply declare them as
%global parameters in both the main script (here) and
%the function itself.  We will treat y as
%a vector with y(1) = y and y(2) = y'.  
pause
global a w
a = 1.5;
w = 2;
dt = 0.1;
%To begin, let's see what the function is actually
%supposed to look like for this set of parameters.
%The exact answer is given in the function yexact.m
%So:
time=[0:.2:20];
v=[0,20,-2,2];
axis(v);
figure(1)
plot(time,yexact(time))
set(gca,'FontSize',18)
pause
hold on
%Now for the Euler method:
y(1)=1;
y(2)=0;
t=0;
i=1;
figure(1)
plot(t,y(1),'o')
set(gca,'FontSize',18)
while t<20,
  y=y+dt*fa(y,t);
  t=t+dt;
  plot(t,y(1),'o')
  set(gca,'FontSize',18)
  erem(i)=y(1)-yexact(t);
  i=i+1;
  echo off
end
echo on
hold off
drawnow
pause
%As you can see, the Euler method starts out following
%the function quite well, but the instability sets
%in at larger times and the error grows quite large.
%The error is kept in the variable erem:
figure(2)
plot(erem)
set(gca,'FontSize',18)
%As you can see, the error gets pretty big!
pause

%We can also look at the Backward Euler method.  In
%this case the method is implicit, so we have to solve
%the root to an equation at each step.  Thus:
figure(1)
plot(time,yexact(time))
axis(v);
set(gca,'FontSize',18)
hold on
y(1)=1;
y(2)=0;
t=0;
i=1;
[m,n] = size(y);
tol=.0001;
eps=.00001;
figure(1)
plot(t,y(1),'o')
set(gca,'FontSize',18)
drawnow
%Because of the modification of the echo command, we will
%have to turn echo off for this routine.  Look at the 
%source code to see how we solve this problem
pause
echo off
while t<20,
  dy=1;
  ybe=y+dt*fa(y,t); %We start with the EM estimate.
  while sum(abs(dy))>tol,
    %This is the function which should be zero:
    g=ybe-y-dt*fa(ybe,t+dt);
    for k=1:n
      e=zeros(1,n);
      e(k)=eps;
      j(1:n,k)=(fa(ybe+e,t+dt)-fa(ybe-e,t+dt))'/2/eps;
    end
    %This is the derivative of the above function:
    gp=eye(n)-dt*j;
    %Now we solve for the updated BE prediction:
    dy=gp\g';
    ybe=ybe-dy';
  end
  t=t+dt;
  y=ybe;
  plot(t,y(1),'o')
  set(gca,'FontSize',18)
  erbe(i)=y(1)-yexact(t);
  i=i+1;
end
echo on
hold off
drawnow
pause
%This tracks the function about the same.  The
%magnitude of the error stays bounded, but still
%grows with time.  We can plot the two errors 
%together:
figure(2)
plot(erem)
set(gca,'FontSize',18)
hold on
plot(erbe,'r')
set(gca,'FontSize',18)
hold off
xlabel('iteration')
ylabel('error')
legend('Euler Method','Backward Euler Method')
drawnow
%As you can see, the Backward Euler method is not a
%great improvement.  This is because the differential
%equation itself is not stiff, but rather is only
%neutrally stable without damping.  -Any- finite
%value of dt will lead to numerical instability!
pause

%OK, now let's look at the trapezoidal rule.  This
%method is also implicit, so again we will use
%Newton's method to solve for each new y.  We get:
figure(1)
plot(time,yexact(time))
axis(v);
set(gca,'FontSize',18)
hold on
y(1)=1;
y(2)=0;
t=0;
i=1;
[m,n] = size(y);
tol=.0001;
eps=.00001;
plot(t,y(1),'o')
set(gca,'FontSize',18)
drawnow
%Again, we turn echo off.  Look at the source code
%to see how it works.
pause
echo off
while t<20,
  dy=1;
  ytr=y+dt*fa(y,t);
  while sum(abs(dy))>tol,
    g=ytr-y-dt*(fa(ytr,t+dt)+fa(y,t))/2;
    for k=1:n
      e=zeros(1,n);
      e(k)=eps;
      j(1:n,k)=(fa(ytr+e,t+dt)-fa(ytr-e,t+dt))'/2/eps;
    end
    gp=eye(n)-dt*j/2;
    dy=gp\g';
    ytr=ytr-dy';
  end
  t=t+dt;
  y=ytr;
  figure(1)
  plot(t,y(1),'o')
  set(gca,'FontSize',18)
  ertr(i)=y(1)-yexact(t);
  i=i+1;
end
echo on
hold off
drawnow
pause
%The trapezoidal rule does a much better job of tracking
%the function, as one would expect.  It's error still
%grows with time, however, and the implicit nature
%of the rule (albeit valuable for stiff equations) is
%fairly slow.  Let's compare this rule to one of the
%simplest of the higher order explicit rules, the
%two stage runge kutta rule:
pause
figure(1)
plot(time,yexact(time))
axis(v);
set(gca,'FontSize',18)
hold on
y(1)=1;
y(2)=0;
t=0;
i=1;
plot(t,y(1),'o')
set(gca,'FontSize',18)
while t<20,
  k1=dt*fa(y,t);
  k2=dt*fa(y+k1,t+dt);
  y=y+(k1+k2)/2;
  t=t+dt;
  figure(1)
  plot(t,y(1),'o')
  set(gca,'FontSize',18)
  er2s(i)=y(1)-yexact(t);
  i=i+1;
  echo off
end
echo on
hold off
drawnow
pause
%Here the method also tracks the function pretty well,
%with a slightly greater error.  We can see this a bit
%better if we plot the two errors simultaneously:
figure(2)
plot(er2s,'r')
set(gca,'FontSize',18)
hold on
plot(ertr,'k')
set(gca,'FontSize',18)
legend('RK-2','TR')
xlabel('iteration')
ylabel('error')
hold off
pause
%From this plot it is apparent that the TR error is
%about half that of the 2-stage runge kutta technique.
%The discrepancy would be -much- larger for a very
%stiff problem in which the explicit method would
%fail, but the implicit method would be stable.
pause
%Now we turn to our final rule, the 4-stage runge
%kutta technique.  We have the simple set of
%equations:
figure(1)
plot(time,yexact(time))
axis(v);
set(gca,'FontSize',18)
hold on
y(1)=1;
y(2)=0;
t=0;
i=1;
plot(t,y(1),'o')
set(gca,'FontSize',18)
while t<20,
  k1=dt*fa(y,t);
  k2=dt*fa(y+k1/2,t+dt/2);
  k3=dt*fa(y+k2/2,t+dt/2);
  k4=dt*fa(y+k3,t+dt);
  y=y+(k1+2*k2+2*k3+k4)/6;
  t=t+dt;
  plot(t,y(1),'o')
  set(gca,'FontSize',18)
  er4s(i)=y(1)-yexact(t);
  i=i+1;
  echo off
end
drawnow
echo on
hold off
pause
%This method also tracks the function very well.  The
%error now is very small, yielding:
figure(2)
plot(er4s)
set(gca,'FontSize',18)
drawnow
%Which is several orders of magnitude smaller than
%the error in the two-stage rule!
pause

hold off
echo off