Trying to get a plot but I get the error below, any way to plot without changing the code but using different syntax?

| | The function is:||

function yprime = model(t,y)
d=5*10^-6; 
xHe=0:1:100; %Important 
xAr=100:-1:0; %Important
pg=(1.5979.*xAr)+(.1624.*xHe); 
V=(pi/6)*d^3; 
As=(pi/4)*d^2; 
Ad=(pi*d^2); 
g=9.81; 
pm=3545; 
CM=0.2165;  
yg=1.66; 
Po=2.76; 
Pe=0.0018; 
Ma=CM*(sqrt((2/(yg-1))*(((Po/Pe)^((yg-1)/yg))-1))); 
mmolg=((40.*xAr)+(4.*xHe))./1000;
R=8.314; 
To=300; 
Te=To/(1+((yg-1)/2)*(Ma^2)); 
vs=sqrt(((yg*R*Te)./(mmolg))); 
vgo=Ma.*vs; 
k=(3.04*(10^-4)).*(vgo.^1.24); 
e=0.035; 
s=5.671*10^-8; 
Kg=(.0179.*xAr)+(.15015.*xHe); 
Cpg=(520.76.*xAr)+(5278.*xHe); 
ug=((224.3*10^-6).*xAr)+((198.6*10^-6).*xHe); 
Tgas=370;
yprime(1,1)=y(2);  
CD=(.28+(((6.*(sqrt(((pg.*d.*(abs((vgo.*exp(-y(1)./k))-y(2))))./ug))))+21)./((pg.*d.*(abs((vgo.*exp(-y(1)./k))-y(2))))./ug))); 
yprime(2,1)=(V.*(pm-pg).*g-(0.5.*pg.*As.*CD).*(abs(y(2)-(vgo.*exp(-y(1)./k)))).*(y(2)-vgo.*exp(-y(1)./k)))./(pm*V); 
vg=vgo.*exp(-y(1)./k); 
vd=y(2); 
Re=((pg.*d.*(abs(vg-vd)))./ug); 
Pr=(ug.*Cpg)./Kg; 
h=(Kg./d).*(2+0.6.*(Re.^.5).*(Pr.^(1/3)));
Cpd=(21.8+.009.*y(3)); 
fr=.043;
yprime(3,1)=fr.*((-h.*Ad*(y(3)-Tgas))-(Ad.*e.*s.*((y(3).^4)-Tgas.^4)))/(V.*pm.*Cpd);

_ | |The script to solve the function is:||_

y0(1,1)=0;
  y0(2,1)=0;
  y0(3,1)=1373;

    xHe=0:1:100; %Important 
    xAr=100:-1:0; %Important
    pg=(1.5979.*xAr)+(.1624.*xHe); 
    g=9.81; 
    pm=3545; 
    CM=0.2165;   
    yg=1.66; 
    Po=2.76;  
    Pe=0.0018; 
    Ma=CM*(sqrt((2/(yg-1))*(((Po/Pe)^((yg-1)/yg))-1))); 
    mmolg=((40.*xAr)+(4.*xHe))./1000; 
    R=8.314; 
    To=300; 
    Te=To/(1+((yg-1)/2)*(Ma.^2)); 
    vs=sqrt(((yg*R*Te)./(mmolg))); 
    vgo=Ma.*vs; 
    k=(3.04*(10^-4)).*(vgo.^1.24); 
    Kg=(.0179.*xAr)+(.15015.*xHe); 
    Cpg=(520.76.*xAr)+(5278.*xHe); 
    ug=((224.3*10^-6).*xAr)+((198.6*10^-6).*xHe); 
    Tgas=370; 
    Pr=(ug.*Cpg)./Kg; 
    fr=.043; 
    tspan1=[0,.005];
    [t1,y1]=ode45('model',tspan1,y0);
    d1=5*10^-6; 
    z1=y1(:,1);  
    vd1=y1(:,2); 
    vg1=vgo.*exp(-z1./k); 
    T1=y1(:,3); 
    Re1=((pg.*d1.*(abs(vg1-vd1)))./ug); 
    Cpd1=21.8+(0.009.*T1);
    h1=(Kg/d1).*(2+0.6.*(Re1.^.5).*Pr.^(1/3)); 
    Tdot1=mean((6*fr*(T1-Tgas).*h1)./(pm*Cpd1.*d1));
    semilogy(xHe,Tdot1,'b*')

_ | The error I get is:|_

>> solve_model
Subscripted assignment dimension mismatch.

Error in model (line 31)
yprime(2,1)=(V.*(pm-pg).*g-(0.5.*pg.*As.*CD).*(abs(y(2)-(vgo.*exp(-y(1)./k)))).*(y(2)-vgo.*exp(-y(1)./k)))./(pm*V);

Error in odearguments (line 87)
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.

Error in ode45 (line 113)
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...

Error in solve_model (line 31)
[t1,y1]=ode45('model',tspan1,y0);

Answer 1

[dpb]

xHe=0:1:100; %Important 
xAr=100:-1:0; %Important
pg=(1.5979.*xAr)+(.1624.*xHe); 
...
yprime(2,1)=(V.*(pm-pg).*g-(0.5.*pg.*As.*CD).*...

"You can't do that!!!" You've defined xHe and xAr as vectors and thereby computed a vector result for the second derivative term yprime(2,:) but as you've written it y(2,1) is a single element of an array. Multiple values can't go into a single location so this is impossible.

Per the documentation of all the ode solvers, the functional form expected is for the routine to expect a single time value, t, and the vector of current y at that time and return a vector of length(y), the values of the derivatives at that one time point. So, even if you could do the assignment or rewrote the function to accept that by returning a 2D array it wouldn't be a functional usable by ode45 (or any of the other solvers).

The cure is to remove the vector expression for the two variables and follow the rules as given in the odefun function. If the idea is that these are time-dependent properties in the system being solved, see Example 3 of the documentation for one way in which to solve for time-dependent parameters using interpolation.

Or, if the intent is to plot the time-dependent solution as a function of the two parameter vectors, then you'll have to wrap the solution to the ODE for each combination of values desired inside a loop and save those results for subsequent plotting.

meshgrid could be of value here along with arrayfun, perhaps to construct such a set of solutions.

What the solution to the dilemma is thus depends on what, specifically, is the end objective but it's not quite "just syntax" (well, actually, it could be said anything in coding reduces to that, but this is a logic issue first before getting to the actual implementation details).

[G]

Thanks, yes the objective is to plot the time-dependent solution as a function of the two parameter vectors (more specifically the xHe variable), but I am unfamiliar with meshgrid and arrayfun.

[dpb]

For only the one variable, then you only need a loop over

for xHe=0:100

inside which you call the ODE solver but use another function that accepts the paramterized variable to introduce it into the solution. See the link on "Parameterizing Functions" under the documentation for ODE45 for examples of how this is done.

[G]

Thanks!