solving differential and non-differential equation

I am trying to solve dydt=2*t + A(t) and A(t)=(y). I would like to feed the new value of A back into the equation of dydt at each time step to get the new value of y. The results i am getting using the code below are not right. Anyone has an idea how to go about this? Thanks.

% Define the time span
tspan = [0 20]; % Define the time span from 0 to 20

% Define the initial condition
y0 = 0; % Initial value of y

% Initialize arrays to store results
t_values = []; % Array to store time values
y_values = []; % Array to store y values
A_values = []; % Array to store A values

% Set initial value of y
y_current = y0;

% Define the loop to iterate over each time step
for t_current = tspan(1):0.5:tspan(2) % Using a step size of 0.5
    % Update A based on current value of y
    A_current = sqrt(y_current);
    
    % Solve for the next value of y using updated A
    [~, y_next] = ode45(@(t, y) odefunc(t, y, A_current), [t_current t_current+0.5], y_current);
    
    % Store current time, y, and A values
    t_values = [t_values; t_current];
    y_values = [y_values; y_current];
    A_values = [A_values; A_current];
    
    % Update y_current for the next iteration
    y_current = y_next(end);
end

% Display the values of y and A at each time step
disp('Time     Value of y     Value of A');
disp([t_values, y_values, A_values]);

% Plot the results
subplot(2,1,1);
plot(t_values, y_values, '-o'); % Plot y(t)
xlabel('Time (t)');
ylabel('y');
title('Solution of dy/dt = 2t + sqrt(y)');

subplot(2,1,2);
plot(t_values, A_values, '-o'); % Plot A(t)
xlabel('Time (t)');
ylabel('A');
title('Solution of A(t) = sqrt(y)');


function dydt = odefunc(t, y, A)
    dydt = 2*t + A;
end

[Torsten]

Why don't you simply set

function dydt = odefunc(t, y)
    dydt = 2*t + sqrt(y);
end

?

[Pc]

Hi @Torsten, I did that and it doesn't seem to create the loop to feed A back. if it did the second time step should give y=1. See the code below

% Define the time span
tspan = [0 20]; % Define the time span from 0 to 20

% Initial conditions
y0 = 0; % Initial value of y

% Solve the system of ODEs
[t, y] = ode45(@odefunc, tspan, y0);

% Calculate m as the square root of y
m = sqrt(y);

% Display the values of t,y and m
disp('Time     Value of y     Value of A');

Time     Value of y     Value of A

disp([t, y, m]);

         0         0         0
    0.5000    0.3477    0.5897
    1.0000    1.5298    1.2368
    1.5000    3.6296    1.9051
    2.0000    6.5632    2.5619
    2.4896   10.1689    3.1889
    2.9792   14.5609    3.8159
    3.4688   19.7396    4.4429
    3.9583   25.7048    5.0700
    4.4583   32.6082    5.7104
    4.9583   40.3318    6.3507
    5.4583   48.8756    6.9911
    5.9583   58.2396    7.6315
    6.4583   68.4238    8.2719
    6.9583   79.4282    8.9122
    7.4583   91.2527    9.5526
    7.9583  103.8975   10.1930
    8.4583  117.3624   10.8334
    8.9583  131.6476   11.4738
    9.4583  146.7529   12.1142
    9.9583  162.6784   12.7545
   10.4583  179.4241   13.3949
   10.9583  196.9900   14.0353
   11.4583  215.3761   14.6757
   11.9583  234.5824   15.3161
   12.4583  254.6089   15.9565
   12.9583  275.4556   16.5969
   13.4583  297.1224   17.2372
   13.9583  319.6095   17.8776
   14.4583  342.9168   18.5180
   14.9583  367.0442   19.1584
   15.4583  391.9918   19.7988
   15.9583  417.7597   20.4392
   16.4583  444.3477   21.0796
   16.9583  471.7559   21.7199
   17.4583  499.9843   22.3603
   17.9583  529.0329   23.0007
   18.4688  559.5326   23.6544
   18.9792  590.8869   24.3082
   19.4896  623.0960   24.9619
   20.0000  656.1598   25.6156

% Plot the results
subplot(2,1,1);
plot(t, y, '-o'); % Plot y(t)
xlabel('Time (t)');
ylabel('y');
title('Solution of dy/dt = 2t');

subplot(2,1,2);
plot(t, m, '-o'); % Plot m(t)
xlabel('Time (t)');
ylabel('m');
title('Solution of m = sqrt(y)');

% Define the ODE function

function dydt = odefunc(t, y)
    dydt = 2*t + sqrt(y);
end

[Sam Chak]

@Peace, What is the purpose of the for-loop sequence?

[Torsten]

if it did the second time step should give y=1.

Why ? And what is the second time step ?

[Pc]

i am trying to get the value of A=sqrt(y) at each time step to be used to solve for the next value of y.

[Torsten]

You mean

y(n+1) = y(n) + (t(n+1)-t(n))*(2*t(n) + sqrt(y(n)))

?

But ode45 does this for you - with a much exacter method.

[Pc]

if you could manually work out the first 3 steps. could you help me confirm that your answers match what is returned by the ODE solver?

[Torsten]

t = 0:0.1:20;
n = numel(t);
y = zeros(n,1);
y(1) = 0;
for i = 2:n
    y(i) = y(i-1)+ (t(i)-t(i-1))*(2*t(i-1) + sqrt(y(i-1)));
end
hold on
plot(t,y)
fun = @(t,y)2*t+sqrt(y);
[T,Y] = ode45(fun,[0 20],0);
plot(T,Y)
hold off
grid on

[Pc]

Hi @Torsten, thanks. I would have expected y and Y to have the same values, but its not the case. Also, Shouldn't it be y(i) = y(i-1)+ (t(i)-t(i-1))*(2*t(i) + sqrt(y(i-1)))?

[Torsten]

I would have expected y and Y to have the same values, but its not the case.

ode45 is more precise than one simple line of code that updates the y-values. But I'd say the two curves are close.

Also, Shouldn't it be y(i) = y(i-1)+ (t(i)-t(i-1))*(2*t(i) + sqrt(y(i-1)))?

No. The explicit Euler method to solve the differential equation dy/dt = f(t,y) evaluates the right-hand side at the old value for t and the old values for y:

y(i) = y(i-1) + (t(i)-t(i-1))*f(t(i-1),y(i-1))

[Pc]

thank you @Torsten for that clarifcation.