Inconsistency of the figure and value on the curve

Hi all. I wrote a code to find the curve for trace and determinant. The problem is that when I zoom and choose one value on the curve and put it into the code, the value obtained for T is different from what is shown on the curve.

For example, here, for c=0.05826 on the curve is -4.06892e-07, while when code runs for this specific value, gives me T =-2.4350e-04.

Could you please let me know what the problem is that this happened and how I can solve it? Thanks in advance for any help

vplc=0.19;
delta=2.5;
tau_max=0.18;
Ktau=0.045;
kc=0.1; 
kh=0.05;
Kbar=0.000015;
kp=0.15;
kb=0.4;
Kpm=0.15;
Ke=7;
vs=0.002;
ks=0.1;
kplc=0.055;
ki=2;
gamma=5.5;
kipr=0.18;

c=0.05826
%c =0:0.000001:0.3; %first I plotted the curves for this interval, zoom and choose a point on the curve and run the code for that specific vale
p=(vplc./ki).*(c.^2./((kplc).^2+c.^2));
h=kh.^4./(c.^4+kh.^4);
Po=(p.^2.*c.^4.*h)./(p.^2.*c.^4.*h.*(1+kb)+kb.*kp.^2.*(kc.^4+c.^4-(c.^4.*kh.^4./(c.^4+kh.^4))));
ct=(vs./(gamma.*kipr.*Po)).*(c.^2./(c.^2+ks.^2))+c.*(1+(1/gamma));
Podenominator=(p.^2.*c.^4.*h.*(1+kb)+kb.*kp.^2.*(kc.^4+c.^4-((c.^4.*kh.^4)./(c.^4+kh.^4))));
Dcnumerator=(8.*c.^7.*kh.^4*vplc.^2)./(ki.^2*(c.^4 + kh.^4).*(c.^2 + kplc.^2).^2) - (4.*c.^9.*kh.^4.*vplc.^2)./(ki.^2.*(c.^4 + kh.^4).*(c.^2 + kplc.^2).^3) - (4.*c.^11.*kh.^4*vplc.^2)./(ki.^2.*(c.^4 + kh.^4).^2.*(c.^2 + kplc.^2).^2);
Dcdenominator=kb.*kp.^2.*(4.*c.^3 - (4*c.^3*kh.^4)./(c.^4 + kh.^4) + (4.*c.^7*kh.^4)./(c.^4 + kh.^4).^2) + (8.*c.^7*kh.^4.*vplc.^2.*(kb + 1))./(ki.^2.*(c.^4 + kh.^4).*(c.^2 + kplc.^2).^2) - (4*c.^9.*kh.^4.*vplc.^2.*(kb + 1))/(ki.^2.*(c.^4 + kh.^4).*(c.^2 + kplc.^2).^3) - (4*c.^11.*kh^4.*vplc^2.*(kb + 1))./(ki^2.*(c.^4 + kh^4).^2.*(c.^2 + kplc.^2).^2);
DcPo=(Dcnumerator.*Podenominator-Dcdenominator.*(p.^2.*c.^4.*h))./(Podenominator.^2);
DhPo=((c.^4.*p.^2.*(Podenominator))-(c.^8.*p.^4.*(kh.^4./(c.^4+kh.^4)).*(1+kb)))./(Podenominator.^2);
Dcf1=kipr.*DcPo.*gamma.*ct-kipr.*DcPo.*(1+gamma).*c-kipr.*Po.*(gamma+1)-((2*vs.*c.*(ks^2))./(c.^2+ks^2).^2);
Dhf1=kipr.*DhPo.*(gamma.*ct-(1+gamma).*c);
Dcf2=(4.*c.^3./tau_max).*((kh.^8./(c.^4+kh.^4).^2))-(4.*c.^3./tau_max).*(kh.^4./(c.^4+kh.^4));
Dhf2=-c.^4./tau_max;


T = Dcf1 + Dhf2
D = Dcf1.*Dhf2-Dhf1.*Dcf2
eig1=(T+sqrt(T.^2-4.*D))./2;
eig2=(T-sqrt(T.^2-4.*D))./2;

figure
plot(c,T,'b')
hold on;
plot(c,D,'r')
hold on
plot(c,T.^2-4.*D,'g')
xlabel('c')
ylabel('T & D')
xlim([0 0.3])
ylim([-0.001 0.003])

Answer 1

[Paul]

Hi @M

The equation for Dcdenominator is not properly vectorized. In one spot it uses /, which should be ./

% was
% Dcdenominator=kb.*kp.^2.*(4.*c.^3 - (4*c.^3*kh.^4)./(c.^4 + kh.^4) + (4.*c.^7*kh.^4)./(c.^4 + kh.^4).^2) + (8.*c.^7*kh.^4.*vplc.^2.*(kb + 1))./(ki.^2.*(c.^4 + kh.^4).*(c.^2 + kplc.^2).^2) - (4*c.^9.*kh.^4.*vplc.^2.*(kb + 1))/(ki.^2.*(c.^4 + kh.^4).*(c.^2 + kplc.^2).^3) - (4*c.^11.*kh^4.*vplc^2.*(kb + 1))./(ki^2.*(c.^4 + kh^4).^2.*(c.^2 + kplc.^2).^2);
%
% should be
% xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx^use./

Dcdenominator=kb.*kp.^2.*(4.*c.^3 - (4*c.^3*kh.^4)./(c.^4 + kh.^4) + (4.*c.^7*kh.^4)./(c.^4 + kh.^4).^2) + (8.*c.^7*kh.^4.*vplc.^2.*(kb + 1))./(ki.^2.*(c.^4 + kh.^4).*(c.^2 + kplc.^2).^2) - (4*c.^9.*kh.^4.*vplc.^2.*(kb + 1))./(ki.^2.*(c.^4 + kh.^4).*(c.^2 + kplc.^2).^3) - (4*c.^11.*kh^4.*vplc^2.*(kb + 1))./(ki^2.*(c.^4 + kh^4).^2.*(c.^2 + kplc.^2).^2);

[M]

Excellent @Paul, it works perfectly... thanks

I have a question, I would appreciate if you can answer it.

With another method that I trust, I found a point where T~0 (when c~0.091), but even though all my calculations (DcPo, DhPo, T and D...)and parameters are correct here (I checked them in Maple), it contradicts with what I got (for example here for c~0,058, T becomes almost zero). you don't have any suggestion how to find the problem?

[Paul]

I'm not following ....