MATLAB: How to get the t=infinity values of an ODE against one of the parameters and then plot those infinity values to a different parameter

MATLABode45

At the moment I'm using the following code to get the t=infinity values against B.*1000./gamma
B = 0:0.01:1;
gamma = 0.1;
sigma = 0:0.01:1;
t_inf_values = zeros(size(B));
for i = 1:numel(B)
beta = B(i);
fprintf('beta = %f\n', beta)
f = @(t,x) [-beta.*x(1)*x(2) + sigma.*x(3); beta.*x(1)*x(2)-gamma.*x(2); gamma.*x(2)-sigma.*x(3)];
[t,xa] = ode45(f, [0 1000], [950 50 0]);
t_inf_values(i) = xa(end,3);
t_inf_values1(i) = xa(end,2);
end
plot((B.*1000./gamma), t_inf_values, 'g')
hold on
plot((B.*1000./gamma), t_inf_values1, 'r')
hold off
Now I need to vary sigma and plot sigma against the infinity values of B.*1000./gamma.

Best Answer

You can create a nested for-loop and then create a surf plot between beta, sigma, and x. Note that consider decreasing the step size for B and Sigma from 0.01 to some higher value. It can take a long time to run
B = 0:0.01:1;
gamma = 0.1;
Sigma = 0:0.01:1;
t_inf_values1 = zeros(numel(B), numel(Sigma));
t_inf_values2 = zeros(numel(B), numel(Sigma));
for i = 1:numel(B)
    for j = 1:numel(Sigma)
        beta = B(i);
        sigma = Sigma(j);
        fprintf('beta = %f,\t sigma = %f\n', beta, sigma) 
        f = @(t,x) [-beta.*x(1)*x(2) + sigma.*x(3); beta.*x(1)*x(2)-gamma.*x(2); gamma.*x(2)-sigma.*x(3)];
        [t,xa] = ode45(f, [0 1000], [950 50 0]);
        t_inf_values1(i, j) = xa(end,3);
        t_inf_values2(i, j) = xa(end,2);
    end
end
subplot(2,1,1);
surf(B, Sigma, t_inf_values1);
shading interp
xlabel('beta');
ylabel('Sigma');
zlabel('x_3(t_{inf})');
subplot(2,1,2);
surf(B, Sigma, t_inf_values2);
shading interp
xlabel('beta');
ylabel('Sigma');
zlabel('x_2(t_{inf})');
Related Question