3 dimensional plotting by varying the initial conditions

I have a set of four differential equations(DE) and I want to see how one of them varies with the other two. I can solve the equation and use plot3 to see the trajectory, but i want to do more and get more information out of it by some kind of 3d plots. I am stuck and hence don't have a lot to show my attempt but here is what i did. The four variables in the system are N, T, A and I. I want to vary the initial conditions of only two of them in some way and solve the DE by fixing the third one. Is it possible i can get a 3d plot of say, how N changes with variation in T and I.

I can change the IC completely and use a for loop and get a set of lines in 3D but i was looking for some sort of 3D surface plot.

My attempt: (Honestly this is not much of an attempt but rather what I have done)

timerange= 0:0.5:200;
IC= [1,0.1,0,0];%initial conditions
[t,y] =ode45(@(t,y) fn(t,y),timerange, IC);

% For 3d plotting defining a meshgrid to specify Initial condition IC
T0=linspace(0,1,100);
I0=linspace(0,1,100);

[T0,I0]= meshgrid(T0,I0);

%Can't figure out how to procceed from here.
plot3(y(1,:),y(3,:),y(2,:));
xlabel('N')
ylabel('I')
zlabel('T')
grid on
xlabel('Time')
ylabel('Fraction of Population')
title('Evolution of cells')
grid on

function rk1 =fn(~,y)

r      = 0.60516;
K      = 8.97523;
A0     = 0.4;
gammaA = 0.04;
eps    = 0.00379;
rho    = 0.02733;
alpha1 = 2.15877;
c1     = 0.02718;
I0     = 0.3;
gammaI = 0.0208;
m      = 0.08006;
q      = 0.61760;gammaT=4;
T0     = 0.4;
alpha  = (gammaA*m)/(A0);
beta   = (gammaT*q)/(T0);


n= y(1);
A= y(2);
T= y(3);
I= y(4);
    

rk1(1)= r*n*(1- n/K)+ alpha*A*n -eps*n*I -beta*n*T;
rk1(2) = A0 - gammaA*A;
rk1(3) = T0- gammaT*T;
rk1(4) = I0 + (rho*I*n)/(alpha1+n) -c1*I*n - gammaI*I;
rk1=rk1(:);
end

Any sort of help or an hint would be very helpful. Thank you.

Answer 1

[darova]

Try quiver3 and streamline

clc,clear
cla
r      = 0.60516;
K      = 8.97523;
A0     = 0.4;
gammaA = 0.04;
eps    = 0.00379;
rho    = 0.02733;
alpha1 = 2.15877;
c1     = 0.02718;
I0     = 0.3;
gammaI = 0.0208;
m      = 0.08006;
q      = 0.61760;gammaT=4;
T0     = 0.4;
alpha  = (gammaA*m)/(A0);
beta   = (gammaT*q)/(T0);

A = 1;

[n,T,I] = meshgrid(0:.2:1);

dn = r*n.*(1- n/K)+ alpha*A.*n -eps*n.*I -beta*n.*T;
% rk1(2) = A0 - gammaA*A;
dT = T0 - gammaT*T;
dI = I0 + (rho*I.*n)./(alpha1+n) -c1*I.*n - gammaI*I;

[X,Z] = meshgrid(0:.2:1);
quiver3(n,T,I,dn,dT,dI,2)
streamline(n,T,I,dn,dT,dI,X,X*0+1,Z)
axis vis3d

xlabel('n axis')
ylabel('T axis')
zlabel('I axis')

[Vira Roy]

Thanks for the reply Darova and I really appreciate the effort. But I suppose the equations that you have used in your modified code is now no longer a differential equations which was governing the system. It is more like an expression now. Since this a sort of tumor model I had to work with the solutions of the differential equations. Your earlier much appreciated help (if you remember) was nice , as you used the differential equation there.

Here, my some sort of aim was to change the initial conditions at finer level and get a 3d plot as a surface to see how the solution of any one of them (say T) changes with other two. It may be a naive question and may be doesn't make sense as well, but if you can help in this, i would be helpful.

Thanks Darova anyways.

[darova]

try this

[Vira Roy]

Thanks for you kind effort Darova, but unfortunately this was not what I was looking for. I had to work with the solutions obtained from differential equations not the expressions . Still I am really thankful that you took time to try the problem.

[darova]

I did my best