%  This function solves homework problem [H7.1]
%  --------------------------------------------
%
echo on
diary nsds_hw7_1.dry
%
%  Set initial conditions
%
t = 0;
tf = 50e-6;
x = [0;0];
y = out_eq2(x,t);
h = 0.5e-6;
tvec = t;
yvec = y;
hvec = h;
tol = 1e-5;
errl = 0;
%
%  Simulate using RKF4/5
%
while t < tf,
  [x4,x5] = rkf45_step2(x,t,h);
  err = norm(x4 - x5,'inf')/max([norm(x4),norm(x5),1.0e-10]);
  if err > tol,
    %
    %  Reject current step
    %
    h = (0.8*tol/err)^(0.2) * h;
    errl = 0;
   else
    %
    %  Accept current step
    %
    t = t + h;
    x = x5;
    y = out_eq2(x,t);
    tvec = [tvec,t];
    yvec = [yvec,y];
    if errl > 0,
      h = (0.8*tol/err)^(0.06) * (errl/err)^(0.08) * h;
     else
      h = (0.8*tol/err)^(0.2) * h;
    end,
    hvec = [hvec,h];
    errl = err;
  end,
end
%
%  Plot the results
%
subplot(3,1,1)
plot(tvec,yvec(1,:))
grid on
title('Homework [H7.1]')
xlabel('Time')
ylabel('u3')
subplot(3,1,2)
plot(tvec,yvec(2,:))
grid on
xlabel('Time')
ylabel('iC')
subplot(3,1,3)
plot(tvec,hvec)
grid on
xlabel('Time')
ylabel('h')
%
print -deps nsds_hw7_1.eps
diary off
echo off
%
return