Matlab Differential Equations Euler’s method

Question

I need help plotting a differential equation ... it keeps coming out all funky and the graph is not what it's supposed to look like.

function [dydt] = diff(y,t)

dydt = (-3*y)+(t*(exp(-3*t)));

end

tI = 0;
yI = -0.1;
tEnd = 5;
dt = 0.5;

t = tI:dt:tEnd;
y = zeros(size(t));
y(1) = yI;

for k = 2:numel(y)
    yPrime = diff(t(k-1),y(k-1));
    y(k) = y(k-1) + dt*yPrime;
end

plot(t,y)
grid on
title('Engr')
xlabel('Time')
ylabel('y(t)')
legend(['dt = ' num2str(dt)])

That's my code, but the graph is not anything like what it's supposed to look like. Am I missing something like an index for the for statement?

Edit

I am getting an error:

Error using diff
Difference order N must be a positive integer scalar.

Error in diff3 (line 12)
    yPrime = diff(t(k-1),y(k-1));

You have problems with function names: matlab already have build-in function 'diff' and try to call this function instead your own function. Try to use different name for your function 'diff'. — Danil Asotsky
Your other problem is that your ODE function is defined as diff(y,t) but you're calling with the arguments flipped: diff(t(k-1),y(k-1)). I'd recommend following what Matlab's ODE solvers do (e.g., ode45) and defining your ODE function's arguments as t, and then y. You'll also find out that your time-step, dt, is too large, even for this simple function. Try 0.1 or less. — horchler
You're also going to want your ODE function to be a sub-function (or a separate M-file function) rather than being that the beginning of your code as you have it shown here. — horchler

Unknown Unknown · Accepted Answer · 2015-01-05T21:13:19

After fixing the errors pointed out by Danil Asotsky and horchler in the comments:

avoiding name conflict with built-in function 'diff'
changing the order of arguments to t,y.
decreasing the time-step dt to 0.1
converting ODE right-hand side to an anonymous function

(and removing unnecessary parentheses in the function definition), your code could look like this:

F = @(t,y) -3*y+t*exp(-3*t);

tI = 0;
yI = -0.1;
tEnd = 5;
dt = 0.1;

t = tI:dt:tEnd;
y = zeros(size(t));
y(1) = yI;

for k = 2:numel(y)
    yPrime = F(t(k-1),y(k-1));
    y(k) = y(k-1) + dt*yPrime;
end

plot(t,y)
grid on
title('Engr')
xlabel('Time')
ylabel('y(t)')
legend(['dt = ' num2str(dt)])

which performs as expected:

output

Matlab Differential Equations Euler’s method

1 Answers