Improved Euler Method with embedded error estimate and variable time step

Greetings all,

I am trying to simulate a system of equations using Improved Euler's method. I'm not really sure if I am doing the embedded error estimate or the variable time step correctly. But I want my code to repeat the previous step through my for loop when e_norm > Emax but am unsure how to do this. Any help is appreciated.

e = 0;
x = .1;
y = .1;
z = .1;

T = .1;
N = 30;

a = 15;

Emax = .1;
Emin = .01;

for n = 1:N+1

      xdot1 = -a*(-0.07*e+0.085*( abs(e+1) - abs(e-1) ));
      ydot1 = -1*(-0.07*e+0.085*( abs(e+1) - abs(e-1) )) - z;
      zdot1 = y;

      xh = x + T*xdot1;
      yh = y + T*ydot1;
      zh = z + T*zdot1;

      xdot2 = -a*(-0.07*e+0.085*( abs(e+1) - abs(e-1) ));
      ydot2 = -1*(-0.07*e+0.085*( abs(e+1) - abs(e-1) )) - zh;
      zdot2 = yh;

      x = x + .5*T*(xdot2 + xdot1);
      y = y + .5*T*(ydot2 + ydot1);
      z = z + .5*T*(zdot2 + zdot1);

      e = y - x;

      X(n) = x;
      Y(n) = y;
      Z(n) = z;
      E(n) = e;

      Error_x(n) = x - xh;
      Error_y(n) = y - yh;
      Error_z(n) = z - zh;

      nx = norm(Error_x,1);
      ny = norm(Error_y,1);
      nz = norm(Error_z,1);

      e_norm = sqrt(nx.^2 + ny.^2 + nz.^2);

      if e_norm > Emax;
          T = T/2;
      elseif e_norm > Emin && e_norm < 0.8*Emax;
          T = 1.05*T;
      elseif e_norm < Emin;
          T = 1.25*T;
      elseif e_norm > 0.8*Emax && e_norm < Emax;
          T = T;
      end

      timestep(n) = T
      Error(n) = e_norm

end

plot(timestep)

% subplot(4,1,1)
% plot(t,X)
% subplot(4,1,2)
% plot(t,Y)
% subplot(4,1,3)
% plot(t,Z)
% subplot(4,1,4)
% plot(t,E)
% 
% figure
% 
% subplot(3,1,1)
% plot(t,Error_x)
% subplot(3,1,2)
% plot(t,Error_y)
% subplot(3,1,3)
% plot(t,Error_z)
% 
figure
% 
% subplot(3,1,1)
%plot(t,nx)
% subplot(3,1,2)
plot(Error)
% subplot(3,1,3)
% plot(t,nz)
% 
% figure
% 
% plot(X,Y)

Answer 1

[Richard Brown]

you could include a while loop inside your inner loop that's loosely like this:

e_norm = inf;
while e_norm > E_max
   do some stuff
   e_norm = something better, hopefully
   if "I've done this too many times and it's not working"
      error('bad')
   end
end