Plotting the function of single variable

I am trying to plot the function N_21(p,s,t) given below.

The code is the following

%--------------------------------------------------------------------------------------------------------------------------------------------

syms h k x p q J

n0 = 400;

r0 = 1/n0;

gamma = 1.26;

r=@(p) r0 * gamma.^p;

N_21 =@(p,s,t) r0 .* n0 .* (2 .* s ./ (2.*s - 1)).^(p+1) .* symprod( r(J), J,[0, p-1]).*...

(exp(- r0 .* (2.*s-1) .* t)./(symprod( r(J) - r0, J,[1,p]))+...

symsum( exp(- r(q) .* (2.*s - 1) .* t)./((r0-r(q)).* (symprod(r(h) - r(q), h, 0, q-1) .* ...

symprod(r(h) - r(q), h, q+1, p))), q,[1, p]));

p2=plot(t, N_21(1,0.85,t),'b--', 'LineWidth',1.5 );

%--------------------------------------------------------------------------------------------------------------------------------------------

"p" and "s" are the parameters and can take positive values (p=1,2,3,4,.... and 0<=s<1). The reason for putting p and s in the code is to be free in choosing any desired value to them.

I can't figure out the error in the code why this is not running. Any suggestion/solution, please.

[John D'Errico]

What error do you get in the code? When you just tell us you got an error, there are many things you may have done wrong. This forces us to run your code, when the error message itself would tell us much about what happened. So when you want to ask about an error, show the COMPLETE text of the error message. EVERYTHING IN RED.

[Iftikhar Ahmed]

@John D'Errico

Below is the error I am getting:

"

Error using symengine

Inconsistent assumptions.

Error in sym/symsum (line 70)

rSym = mupadmex('symobj::map',fsym.s,'symobj::symsum',x.s,a.s,b.s);

Error in

Plots_slopes>@(p,s,t)r0.*n0.*(2.*s./(2.*s-1)).^(p+1).*symprod(r(q),q,0,p-1).*(exp(-r0.*(2.*s-1).*t)./(symprod(r(q+0.0000001)-r0,q,0,p))+symsum(exp(-r(q).*(2.*s-1).*t)./((r0-r(q+0.0000001)).*(symprod(r(h)-r(q),h,1,q-1).*symprod(r(h)-r(q),h,q+1,p))),q,0,p))

Error in Plots_slopes (line 62)

p2=plot(t,N_2(1,0.85,t),'r--', 'LineWidth',1.5 );"

Answer 1

[David Hill]

n0 = 400;
r0 = 1/n0;
gamma = 1.26;
p=20;%assumed
s=.7;%assumed
P=1:(p+1);
r=r0*gamma.^(0:p);
N=@(t)r0*n0*(2*s/(2*s-1))^(p+1)*prod(r(1:end-1))*(exp(-r0*(2*s-1)*t)/prod(r(2:end)-r(1))+sum(exp(-r(2:end)*(2*s-1)*t)./(r(1)-r(2:end)).*arrayfun(@(x)prod(r(P~=x)-r(x)),2:(p+1))));
t=0:.1:10;
n=zeros(size(t));
for k=1:length(t)
    n(k)=N(t(k));
end
plot(t,n);