How to replace ode45 with Eulers method?

This code currently works perfectly but I'm being asked to solve the differential equation using Euler's method instead of ode45 and I'm not sure how to do that. Could anyone help me replace ode45 with Euler's method? Thanks!

clear 
  close all
  clc

    tic

    % INPUT SECTION ========================================================

    % External current stimulus  [Iext = 1.0]
    Iext = 0.33;
    % Initial conditions: y0 = [ v(0) = -2.8, w(0) = -1.8] 
    y0 = [-2.8; -1.8];
    % Simulation time [tMin = 0 tMax = 200]
    tSpan = [0 200];
    % Phase space plot dimensions: X axis [-3 3] / Y axis [-2 3]
    Lx = [-3 3]; Ly = [-2 3];
    % Phase plot ticks
    ticksX = Lx(1):1:Lx(2);
    ticksY = Ly(1):1:Ly(2); 
    % Phase space setup for vector field: number of vectors nX [20] 
    nX = 20;
    % Fitzhugh-Nagoma model parameters
    a = 0.7; b = 0.8; c = 12.5;

    % CALCULATION SECTION ===================================================
    % K use to pass variables to function FNode
       K = [Iext; a; b; c];

    % Solve differential equations using ode45
      [t,y] = ode45(@(t,y) FNode(t,y,K), tSpan,y0);

    % Vector field 
      x1 = linspace(Lx(1),Lx(2),nX);
      x2 = linspace(Ly(1),Ly(2),nX);
      [xx, yy] = meshgrid(x1,x2);
      f = (xx - xx.^3/3 - yy + Iext);   
      g = (1/c).*(xx + a - b.*yy); 
      fs = f./sqrt(f.^2 + g.^2);    % unit vectors
      gs = g./sqrt(f.^2  +g.^2);

    % Critical point - output to Command Window
    syms p
    Sp = vpasolve(p-p^3/3-(p+a)/b + Iext == 0,p,[-3 3]);
    Sq = (Sp+a)/b;
    Sp = double(Sp); Sq = double(Sq);
    disp('Critical point');
    fprintf('   v_C =  %2.2f\n', Sp);
    disp('   ')
    fprintf('   w_C =  %2.2f\n', Sq);

    % GRAPHICS SECTION=======================================================  

       FS = 14;  % fontsize

    % Phase space plot: v vs w   ------------------------------------------
    figure(1)
       pos = [0.35 0.05 0.29 0.39];
       set(gcf,'Units','normalized');
       set(gcf,'Position',pos);
       set(gcf,'color','w');
       hold on
       box on

    % VECTOR FIELD
       hq = quiver(xx,yy,fs,gs);
       set(hq,'color',[0.2 0.2 0.2],'AutoScaleFactor',0.6);
       set(gca,'fontsize',FS)
       xlim(Lx)
       ylim(Ly)
       set(gca,'xtick',ticksX);
       set(gca,'ytick',ticksY);
       grid on   
       xlabel('membrane potential v'); ylabel('recovery variable  w');

    % v nullcline
         v = linspace(Lx(1),Lx(2),200);
         w = (v - v.^3/3 + Iext);
            xP = v; yP = w;
              plot(xP,yP,'r','linewidth',1.5)

              % r nullcline 
        w = (v + a)/b;
            xP = v; yP = w;
              plot(xP,yP,'m','linewidth',1.5)

    % Phase portrait trajectory  
        xP = y(:,1); yP = y(:,2);
          plot(xP,yP,'b','linewidth',2)
        xP = y(1,1); yP = y(1,2);    % initial conditions: start of trajectory
          Hplot = plot(xP,yP,'o');
          set(Hplot,'markersize',8,'markerfacecolor',[0 1 0],'markeredgecolor',[0 1 0])
        xP = y(end,1); yP = y(end,2);   % end of trajectory
          Hplot = plot(xP,yP,'o');
          set(Hplot,'markersize',8,'markerfacecolor',[1 0 0],'markeredgecolor',[1 0 0])

       tm1 = 'I_{ext} = ';
       tm2 = num2str(Iext,'%3.3f');
       tm3 = '  v_C = ';
       tm4 = num2str(Sp,'%3.2f');
       tm5 = '  w_C = ';
       tm6 = num2str(Sq,'%3.2f');
       tm = [tm1 tm2 tm3 tm4 tm5 tm6];   
       hT = title(tm,'FontName','Courier');  

    % Time evolution of v and w   
    figure(2)
      pos = [0.05 0.05 0.29 0.29];
      set(gcf,'Units','normalized');
      set(gcf,'Position',pos);
      set(gcf,'color','w');
      xP = t; yP = y(:,1);  % v
        plot(xP,yP,'b','linewidth',2)
      hold on
        yP = y(:,2);        % w
      plot(xP,yP,'r','linewidth',2)

      legend('v','w','location','south','orientation','horizontal')
      xlabel('t')
      ylabel('v & w')
      title(tm,'fontName','Courier')
      grid on
      set(gca,'fontsize',FS)
      box on 

      disp('  ')
      toc  

    % FUNCTIONS ===========================================================

    function dydt = FNode(t,y,K)
       a = K(2); b = K(3); c = K(4); Iext = K(1);
       %Iext = 0.003 * t;
       %if t > 50; Iext = 0.5; end
       %if t > 55; Iext =   0; end 

       dydt = [(y(1) - y(1)^3/3 - y(2) + Iext); (1/c)*(y(1) + a - b*y(2))];

    end

Answer 1

[Torsten]

...
% Solve differential equations using ode45
  [t,y] = euler(@FNode,tSpan,y0,K);
...

function [tsol, ysol] = euler(fun, tspan, y0, K)
  n = floor(20*tspan(2));
  %Compute h from the time grid.
  h = (tspan(2) - tspan(1))/n;
  %How many equations?
  m = length(y0);
  %Initilize t and y
  t = tspan(1);
  y = reshape(y0, 1, m);
  % Store initial conditions
  tsol = t;
  ysol = y;
  % Euler loop
  for i = 1: n
    % Take the Euler step into the temporary variable
    y1 = y + h * fun(t, y, K).';
    % Update the temporary variables
    t = t + h;
    y = y1;
    % Update the solution vector
    tsol = [tsol ; t];
    ysol = [ysol ; y];
  end 
end

Answer 2

[Jan]

I recommend to ask an internet search engine for "Matlab euler method". You will find e.g.:

• https://www.mathworks.com/matlabcentral/answers/278300-matlab-code-help-on-euler-s-method
• https://www.mathworks.com/matlabcentral/answers/130653-code-of-euler-s-method
• http://people.math.sfu.ca/~ralfw/math467w03/matlab/euler_matlab.pdf
• http://tutorial45.com/euler-method-matlab-code/

Does this help already? Try to implement it. It is not tricky to expand this to a 2D y. Post what you have tried so far and ask a specific question in case of problems.

Answer 3

[KSSV]

Read about ode23

[Westin Messer]

Thanks for the tip! Unfortunately, I know about ode23 and that is not Euler's method. Sometimes ode solvers like ode23 and ode45 make hidden assumptions when calculating that you don't know about so I need to use Euler's method to clearly see the iterative loop and how the ode is being solved.

[Torsten]

Explicit or implicit Euler ?

[Westin Messer]

Explicit Euler

[Torsten]

http://tutorial45.com/euler-method-matlab-code/

[Westin Messer]

I've done Euler's method from scratch before. I'm confused as to how to modify the code that was written for ode45. My question isn't how to do Euler's method, it's how to replace ode45 with Euler's method in a script written for ode45, if that makes sense.

[Torsten]

Replace the line

[t,y] = ode45(@(t,y) FNode(t,y,K), tSpan,y0);

by

[t,y] = euler(@FNode,tSpan,y0);

and add

function [tsol, ysol] = euler(fun,tspan,y0)
...

[Westin Messer]

Thanks! Sorry for my ignorance but what goes inside the function?

[Torsten]

Here is what the function should contain:

https://de.mathworks.com/matlabcentral/answers/366717-implementing-forward-euler-method

[Westin Messer]

I'm having trouble implementing that function because it's for a completely different problem.

[Torsten]

No, that's wrong. The function "euler_method" is written for an arbitrary system of ODEs, thus can be immediately be applied for your case.

[Westin Messer]

Okay, where do I define my differential equations then?

[Torsten]

You can keep your function "FNode" as it is.

That's why I wrote

Replace the line

[t,y] = ode45(@(t,y) FNode(t,y,K), tSpan,y0);

by

[t,y] = euler(@FNode,tSpan,y0);

[Torsten]

Ok, you'll additionally have to provide K to "euler":

[t,y] = euler(@FNode,tSpan,y0,K);

[Westin Messer]

I see, the function in the link you sent is:

function [tgrid, Y] = euler_method(fun, y_0, n, T)

But the one you told me to use is:

function [tsol, ysol] = euler(fun,tspan,y0)

they're different, you went from four variables to three. Which one do I need to delete?

[Jan]

@Westin Messer: You can define the time either by a vector, or by a range and a number of steps. This is fairly equivalent.

[Westin Messer]

This is what I have now but I'm getting a lot of errors. Am I getting close?

clear 
close all
clc

tic

% INPUT SECTION ========================================================

% External current stimulus  [Iext = 1.0]
Iext = 0.33;
% Initial conditions: y0 = [ v(0) = -2.8, w(0) = -1.8] 
y0 = [-2.8; -1.8];
% Simulation time [tMin = 0 tMax = 200]
tSpan = [0 200];
% Phase space plot dimensions: X axis [-3 3] / Y axis [-2 3]
Lx = [-3 3]; Ly = [-2 3];
% Phase plot ticks
ticksX = Lx(1):1:Lx(2);
ticksY = Ly(1):1:Ly(2); 
% Phase space setup for vector field: number of vectors nX [20] 
nX = 20;
% Fitzhugh-Nagoma model parameters
a = 0.7; b = 0.8; c = 12.5;

% CALCULATION SECTION ===================================================
% K use to pass variables to function FNode
   K = [Iext; a; b; c];

% Solve differential equations using ode45
  [t,y] = euler(@FNode,tSpan,y0);

% Vector field 
  x1 = linspace(Lx(1),Lx(2),nX);
  x2 = linspace(Ly(1),Ly(2),nX);
  [xx, yy] = meshgrid(x1,x2);
  f = (xx - xx.^3/3 - yy + Iext);   
  g = (1/c).*(xx + a - b.*yy); 
  fs = f./sqrt(f.^2 + g.^2);    % unit vectors
  gs = g./sqrt(f.^2  +g.^2);

% Critical point - output to Command Window
syms p
Sp = vpasolve(p-p^3/3-(p+a)/b + Iext == 0,p,[-3 3]);
Sq = (Sp+a)/b;
Sp = double(Sp); Sq = double(Sq);
disp('Critical point');
fprintf('   v_C =  %2.2f\n', Sp);
disp('   ')
fprintf('   w_C =  %2.2f\n', Sq);

% GRAPHICS SECTION=======================================================  

     FS = 14;  % fontsize

% Phase space plot: v vs w   ------------------------------------------
figure(1)
   pos = [0.35 0.05 0.29 0.39];
   set(gcf,'Units','normalized');
   set(gcf,'Position',pos);
   set(gcf,'color','w');
   hold on
   box on

% VECTOR FIELD
   hq = quiver(xx,yy,fs,gs);
   set(hq,'color',[0.2 0.2 0.2],'AutoScaleFactor',0.6);
   set(gca,'fontsize',FS)
   xlim(Lx)
   ylim(Ly)
   set(gca,'xtick',ticksX);
   set(gca,'ytick',ticksY);
   grid on   
   xlabel('membrane potential v'); ylabel('recovery variable  w');

% v nullcline
     v = linspace(Lx(1),Lx(2),200);
     w = (v - v.^3/3 + Iext);
        xP = v; yP = w;
          plot(xP,yP,'r','linewidth',1.5)

            % r nullcline 
      w = (v + a)/b;
          xP = v; yP = w;
            plot(xP,yP,'m','linewidth',1.5)

% Phase portrait trajectory  
    xP = y(:,1); yP = y(:,2);
      plot(xP,yP,'b','linewidth',2)
    xP = y(1,1); yP = y(1,2);    % initial conditions: start of trajectory
      Hplot = plot(xP,yP,'o');
      set(Hplot,'markersize',8,'markerfacecolor',[0 1 0],'markeredgecolor',[0 1 0])
    xP = y(end,1); yP = y(end,2);   % end of trajectory
      Hplot = plot(xP,yP,'o');
      set(Hplot,'markersize',8,'markerfacecolor',[1 0 0],'markeredgecolor',[1 0 0])

     tm1 = 'I_{ext} = ';
     tm2 = num2str(Iext,'%3.3f');
     tm3 = '  v_C = ';
     tm4 = num2str(Sp,'%3.2f');
     tm5 = '  w_C = ';
     tm6 = num2str(Sq,'%3.2f');
     tm = [tm1 tm2 tm3 tm4 tm5 tm6];   
     hT = title(tm,'FontName','Courier');  

% Time evolution of v and w   
figure(2)
  pos = [0.05 0.05 0.29 0.29];
  set(gcf,'Units','normalized');
  set(gcf,'Position',pos);
  set(gcf,'color','w');
  xP = t; yP = y(:,1);  % v
    plot(xP,yP,'b','linewidth',2)
  hold on
    yP = y(:,2);        % w
  plot(xP,yP,'r','linewidth',2)

    legend('v','w','location','south','orientation','horizontal')
    xlabel('t')
    ylabel('v & w')
    title(tm,'fontName','Courier')
    grid on
    set(gca,'fontsize',FS)
    box on 

    disp('  ')
    toc  

% FUNCTIONS ===========================================================

function dydt = FNode(t,y,K)
   a = K(2); b = K(3); c = K(4); Iext = K(1);
   %Iext = 0.003 * t;
   %if t > 50; Iext = 0.5; end
   %if t > 55; Iext =   0; end 

     dydt = [(y(1) - y(1)^3/3 - y(2) + Iext); (1/c)*(y(1) + a - b*y(2))];

end

function [tsol, ysol] = euler(fun, tspan, y0)
if nargin(fun) ~=2
    error('fun must take two inputs, t and y.');
end 
if ~all(size(y0) == size(fun(0, y0)))
    error('You have not passed appropriate fun or y_0.');
end 
%Set up the time grid. ***NOTE THE n+1***
tsol = linspace(0, T, n+1);
%Compute h from the time grid.
h = tsol(2) - tsol(1);
%Orient tgrid as a column vector.
tsol = reshape(tsol, n+1, 1);
%How many equations?
m = length(y0);
%Orient y0 as a row vector.
y0 = reshape(y0, 1, m);
% Preallocate an array to hold the approximate solution. Each row
% corresponds to a point in the time grid.
ysol = zeros(n+1, m);
% Set the initial conditions.
ysol(1,:) = y0;
% Euler loop
for i = 1: n
      % Store the point in time as a temporary variable
      t_i = tsol(i);
      % Take the Euler step into the temporary variable
      y_1 = y0 + h * fun(t_i, y0);
      % Store the Euler step
      ysol(i+1,:) = y_1;
      % Update the temporary variable
      y0 = y_1;
end
end