Implementing Euler's Method in GNU Octave

Question

I am reading "Numerical Methods for Engineers" by Chapra and Canale. In it, they've provided pseudocode for the implementation of Euler's method (for solving ordinary differential equations). Here is the pseucode:

Pseucode for implementing Euler's method

I tried implementing this code in GNU Octave, but depending on the input values, I am getting one of two errors:

1. The program doesn't give any output at all. I have to press 'Ctrl + C' in order to break execution.
2. The program gives this message:

error: 'ynew' undefined near line 5 column 21
error: called from
    Integrator at line 5 column 9
    main at line 18 column 7

I would be very grateful if you could get this program to work for me. I am actually an amateur in GNU Octave. Thank you.

Edit 1: Here is my code. For main.m:

%prompt user
y = input('Initial value of y:');
xi = input('Initial value of x:');
xf = input('Final value of x:');
dx = input('Step size:');
xout = input('Output interval:');

x = xi;
m = 0;
xpm = x;
ypm = y;

while(1)
    xend = x + xout;
    if xend > xf
      xend = xf;
      h = dx;
      Integrator(x,y,h,xend);
      m = m + 1;
      xpm = x;
      ypm = y;
      if x >= xf
        break;
      endif
    endif  
end

For Integrator.m:

function Integrator(x,y,h,xend)
    while(1)
      if xend - x < h
        h = xend - x;
        Euler(x,y,h,ynew);
        y = ynew;
        if x >= xend
          break;
        endif
      endif    
    end
endfunction

For Euler.m:

function Euler(x,y,h,ynew)
   Derivs(x,y,dydx);
   ynew = y + dydx * h;
   x = x + h;
endfunction

For Derivs.m:

function Derivs(x,y,dydx)
    dydx = -2 * x^3 + 12 * x^2 - 20 * x + 8.5;
endfunction

Edit 2: I shoud mention that the differential equation which Chapra and Canale have given as an example is:

y'(x) = -2 * x^3 + 12 * x^2 - 20 * x + 8.5

That is why the 'Derivs.m' script shows dydx to be this particular polynomial.

Answer 1

Your functions should look like

  function [x, y] = Integrator(x,y,h,xend)
    while x < xend
      h = min(h, xend-x)
      [x,y] = Euler(x,y,h);
    end%while
  end%function

as an example. Depending on what you want to do with the result, your main loop might need to collect all the results from the single steps. One variant for that is

  m = 1;
  xp = [x];
  yp = [y];  
  while x < xf
    [x,y] = Integrator(x,y,dx,min(xf, x+xout));
    m = m+1;
    xp(m) = x;
    yp(m) = y;
  end%while

Answer 2

Here is my final code. It has four different M-files:

1. main.m

%prompt the user
y = input('Initial value of y:');
x = input('Initial value of x:');
xf = input('Final value of x:');
dx = input('Step size:');
xout = dx;

%boring calculations
m = 1;
xp = [x];
yp = [y];  
while x < xf
  [x,y] = Integrator(x,y,dx,min(xf, x+xout));
  m = m+1;
  xp(m) = x;
  yp(m) = y;
end

%plot the final result
plot(xp,yp);
title('Solution using Euler Method');
ylabel('Dependent variable (y)');
xlabel('Independent variable (x)');
grid on;

2. Integrator.m

%This function takes in 4 inputs (x,y,h,xend) and returns 2 outputs [x,y]
function [x,y] = Integrator(x,y,h,xend)
  while x < xend
    h = min(h, xend-x);
    [x,y] = Euler(x,y,h);
  end
endfunction  

3. Euler.m

%This function takes in 3 inputs (x,y,h) and returns 2 outputs [x,ynew] 
function [x,ynew] = Euler(x,y,h)
   dydx = Derivs(x,y);
   ynew = y + dydx * h;
   x = x + h;
endfunction

4. Derivs.m

%This function takes in 2 inputs (x,y) and returns 1 output [dydx]
function [dydx] = Derivs(x,y)
    dydx = -2 * x^3 + 12 * x^2 - 20 * x + 8.5;
endfunction