Plotting intermediate variables in ode45 solver

Hello, I have a system of three differential equations that I define in the function differential:

function dvdt = differential(t,v)
global rho kappa phi tau xi N1 N2 p b delta alpha mu sigma c_t theta lambda
dvdt = zeros(3,1);
g=v(1); x1=v(2); x2=v(3);
dvdt(1) = rho - kappa*(N1*(1-x1)*g1 + N2*(1-x2)*g2);
dvdt(2) = phi*x1*(pi_t1-pi_1);
dvdt(3) = phi*x2*(pi_t2-pi_2);

The intermediate variables that are used for solving the above three differential equations are as below:

g1 = min(g/(kappa*(N1*(1-x1)+N2*(1-x2))),gstar);
g2 = min(g/(kappa*(N1*(1-x1)+N2*(1-x2))),gstar);
c_g = max(tau*(xi-g),0) + delta;
gstar = (p-c_g)/(2*b);
M = alpha*(N1*x1+N2*x2) + mu*(N1*(1-x1)+N2*(1-x2));
R= (1 - theta/(1+(lambda*M))) * (sigma);
tstar = (p-c_t)/(2*b);
t1 = min(R/(N1*(1-x1)+N2*(1-x2)),tstar);
t2 = max((R-N1*x1*t1)/(N2*(1-x2)), 0);
pi_t1 = max((p*t1(:,1) - b*t1.^2 - c_t*t1(:,1) - alpha),0);
pi_t2 = max((p*t2(:,1) - b*t2.^2 - c_t*t2(:,1) - alpha),0);
pi_g1 = max((p*g1(:,1) - b*g1.^2 - c_g*g1(:,1) - mu),0);
pi_g2 = max((p*g2(:,1) - b*g2.^2 - c_g*g2(:,1) - mu),0);
pi_1=pi_t1*x1+pi_g1*(1-x1);
pi_2=pi_t2*x2+pi_g2*(1-x2);

Currently, I call the differential function in main function and plot the differential functions in the main function:

[T, V] = ode45(@(t,v) differential(t,v), [0:1/200:20], [G0 X1_init X2_init]);
gt = V(:, 1); x1t = V(:, 2);x2t = V(:, 3); 
plot(x1t, x2t, 'k');

My problem is that I need to plot a few of the intermediate variables too. So, I am confused about where to define them, i.e., whether in the main function or in the differential function and how to store/call these variables. I would ideally like to have the code for plotting both the differential functions as well as the intermediate variables in the main function. But, I cannot understand if that is possible. Can anyone help advise? Thank you!

[Sam K]

Any help please??

Answer 1

[Walter Roberson]

See http://www.mathworks.com/help/matlab/math/parameterizing-functions.html

You could pass in a parameter that is a function handle to a function you would call within differential, passing in the values to be plotted. For example you could do

   apfig = figure();
   apax = axes('Parent', pfig);
   ALh = [animatedline('Parent', apax, 'DisplayName', 'pi_t1'), animatedline('Parent', apax, 'DisplayName', 'pi_t2'), animatedline('Parent', apax, 'DisplayName', 'pi_g1'), animatedline('Parent', apax, 'DisplayName', 'pi_g2')];
   auxplot = @(X,Y) arrayfun(@(L,t,x) addpoints(L,t,x), ALh, X, Y);

   [T, V] = ode45(@(t,v) differential(t, v, auxplot, rho, kappa, phi, tau, xi, N1, N2, p, b, delta, alpha, mu, sigma, c_t, theta, lambda), [0:1/200:20], [G0 X1_init X2_init]);

gt = V(:, 1); x1t = V(:, 2);x2t = V(:, 3); 

mpfig = figure();
mpax = axes('Parent', mpfig);
plot(mpax, x1t, x2t, 'k');

with

function dvdt = differential(t, v, auxplot, rho, kappa, phi, tau, xi, N1, N2, p, b, delta, alpha, mu, sigma, c_t, theta, lambda)

dvdt = zeros(3,1);

g1 = min(g/(kappa*(N1*(1-x1)+N2*(1-x2))),gstar);
g2 = min(g/(kappa*(N1*(1-x1)+N2*(1-x2))),gstar);
c_g = max(tau*(xi-g),0) + delta;
gstar = (p-c_g)/(2*b);
M = alpha*(N1*x1+N2*x2) + mu*(N1*(1-x1)+N2*(1-x2));
R= (1 - theta/(1+(lambda*M))) * (sigma);
tstar = (p-c_t)/(2*b);
t1 = min(R/(N1*(1-x1)+N2*(1-x2)),tstar);
t2 = max((R-N1*x1*t1)/(N2*(1-x2)), 0);
pi_t1 = max((p*t1(:,1) - b*t1.^2 - c_t*t1(:,1) - alpha),0);
pi_t2 = max((p*t2(:,1) - b*t2.^2 - c_t*t2(:,1) - alpha),0);
pi_g1 = max((p*g1(:,1) - b*g1.^2 - c_g*g1(:,1) - mu),0);
pi_g2 = max((p*g2(:,1) - b*g2.^2 - c_g*g2(:,1) - mu),0);
pi_1=pi_t1*x1+pi_g1*(1-x1);
pi_2=pi_t2*x2+pi_g2*(1-x2);

auxplot([t t t t], [pi_t1, pi_t2, pi_g1, pi_g2]);

g=v(1); x1=v(2); x2=v(3);
dvdt(1) = rho - kappa*(N1*(1-x1)*g1 + N2*(1-x2)*g2);
dvdt(2) = phi*x1*(pi_t1-pi_1);
dvdt(3) = phi*x2*(pi_t2-pi_2);