Trying to make euler's method program work

1 view (last 30 days)
ce on 5 Dec 2022
Edited: Torsten on 6 Dec 2022
function [x,y] = euler_for(x0,y0,x_end,h,fcn)
%Solve the ivp
% y' = f(x,y), x0 <= x <= b, y(x0) = y0
%Use euler's method with a stepsize h. The
%user must supply a program with some name, say, deriv,
%and a first line of the form
% function ans = deriv(x,y)
%A sample call would be
% [t,z] = euler_for(t0,z0,b,delta,'deriv')
%
%output:
%the routine eulercls will return two vectors, x and y.
%The vector x will contain the node points
% x(1) = x0, x(j) = x0 + (j-1)*h, j=1,2,...,N
%with
% x(N) <= x_end - h, x(N) + h > x_end - h
%The vector y will contain the estimates of the solution Y
%at the node points in x.
%
n = fix((x_end-x0)/h)+1;
x = linspace(x0,x_end,n);
y = zeros(n,1);
y(1) = y0;
for i = 2:n
y(i) = y(i-1)+h*feval(fcn,x(i-1),y(i-1));
end
I'm not sure what is is meant by the function that I need to supply. I know what function I'm supposed to use, but I don't know how to add it to the program. I tried making a function called fcn but when I run the program it says there is an error in fcn, although it doesn't say what it is. This is my fcn code:
function f = fcn(x,y)
f = (cos(atan(x)))^2;
end

Alan Stevens on 5 Dec 2022
Like this:
% Arbitrary values here - insert your own
x0 = 0;
x_end = 10;
h = 0.1;
y0 = 1;
n = fix((x_end-x0)/h)+1;
x = linspace(x0,x_end,n);
y = zeros(n,1);
y(1) = y0;
for i = 2:n
y(i) = y(i-1)+h*fcn(x(i-1));
end
plot(x,y),grid
xlabel('x'),ylabel('y')
function f = fcn(x) % Your function doesn't depend on y
f = (cos(atan(x)))^2;
end
Torsten on 6 Dec 2022
Edited: Torsten on 6 Dec 2022
x0 = 0;
x_end = 10; % Integration interval [x0,x_end]
h = 0.1; % Step size of integration
fcn = @(x,y) (cos(y))^2; % Y' = cos(Y)^2
y0 = 0; % Y(x0) = 0
[x,y] = euler_for(x0,y0,x_end,h,fcn);
hold on
plot(x,y),grid
xlabel('x'),ylabel('y')
syms y(x)
Ysol(x) = dsolve(diff(y,x)==cos(y)^2,y(0)==0);
fplot(Ysol,[x0 x_end])
hold off
function [x,y] = euler_for(x0,y0,x_end,h,fcn)
n = fix((x_end-x0)/h)+1;
x = linspace(x0,x_end,n);
y = zeros(n,1);
y(1) = y0;
for i = 2:n
y(i) = y(i-1)+h*fcn(x(i-1),y(i-1));
end
end