how can i make graph between real(c) and k for diffferent values of h

tic
format shortG
syms c k 
c11=106.8;
c33=54.57;
c13=9.68;
c15=0.28;
c35=-1.69;
c55=25.05;
r1=2727;
l=2.46;
u=5.66;
r2=7800;
n=2.54;
i=sqrt(-1);
c1=c^2;
a01=c33*c55;
a02=c35^2;
a0=a01-a02;
a11=-2*c15*c33;
a12=2*c13*c35;
a1=(a11+a12)*i;
a111=c33+c55;
a211=(a111*r1);
a212=2*c15*c35;
a213=c11*c33;
a214=c13^2;
a215=2*c13*c55;
a21=vpa(2.1712e+05)*c1;
a22=-a212-a213+a214+a215;
a2=vpa(a21-5248.5);
a31=-2*(c15+c35)*r1;
a32=2*c11*c35;
a33=-2*c13*c15;
a3=vpa(((7690.1)*c1-360.98-5.4208)*i);
a41=r1^2;
a42=-(c55+c11)*r1;
a43=c11*c55;
a44=-c15^2;
a4=vpa(7436529*c1^2-(3.5955e+05)*c1+ 2675.3-0.0784);
S=[a0 a1 a2 a3 a4];
R=roots(S);
s11=R(1);
s1=vpa(s11);
s12=R(2);
s2=vpa(s12);
s13=R(3);
s3=vpa(s13);
s14=R(4);
s4=vpa(s14);
m11=c55*(s1^2)-2*i*c15*s1+r1*(c^2)-c11;
m12=c35*(s1^2)-i*(c13+c55)*s1-c15;
m1=-(m11/m12);
m21=c55*(s2^2)-2*i*c15*s2+r1*(c^2)-c11;
m22=c35*(s2^2)-i*(c13+c55)*s2-c15;
m2=-(m21/m22);
m31=c55*(s3^2)-2*i*c15*s3+r1*(c^2)-c11;
m32=c35*(s3^2)-i*(c13+c55)*s3-c15;
m3=-(m31/m32);
m41=c55*(s4^2)-2*i*c15*s4+r1*(c^2)-c11;
m42=c35*(s4^2)-i*(c13+c55)*s4-c15;
m4=-(m41/m42);
alpha1=(l+2*u)/r2;
alpha=vpa(alpha1,4);
beta1=u/r2;
beta=vpa(beta1,4);
d11=(n^2)*(beta/alpha);
d1=vpa(d11,4);
d21=(n^2)*((1-(beta/alpha))^2)+n*(((c^2)/alpha)-n)+n*(beta/alpha)*(((c^2)/beta)-n);
d2=vpa(d21,4);
d31=(beta/alpha)*(((c^2)/alpha)-n)*(((c^2)/beta)-n);
d3=vpa(d31,4);
P=[d1 0 d2 0 d3];
R1=roots(P);
p11=R1(1);
p1=vpa(p11,4);
p12=R1(2);
p2=vpa(p12,4);
n11=((c^2)/alpha)-n+n*(beta/alpha)*(p1^2);
n12=i*n*(1-(beta/alpha))*p1;
n1=n11/n12;
n21=((c^2)/alpha)-n+n*(beta/alpha)*(p2^2);
n22=i*n*(1-(beta/alpha))*p2;
n2=n21/n22;
k1=i*c15-m1*s1*c35+(i*m1-s1)*c55;
k2=i*c15-m2*s2*c35+(i*m2-s2)*c55;
k3=i*c15-m3*s3*c35+(i*m3-s3)*c55;
k4=i*c15-m4*s4*c35+(i*m4-s4)*c55;
k5=-n*u*(i*n1-p1);
k6=-n*u*(i*n2-p2);
k7=i*c13-m1*s1*c33+(i*m1-s1)*c35;
k8=i*c13-m2*s2*c33+(i*m2-s2)*c35;
k9=i*c13-m3*s3*c33+(i*m3-s3)*c35;
k10=i*c13-m4*s4*c33+(i*m4-s4)*c35;
k11=-n*(i*l-(l+2*u)*p1*n1);
k12=-n*(i*l-(l+2*u)*p2*n2);
for h=2.4
    e1=exp(k*s1*h);
    e2=exp(k*s2*h);
    e3=exp(k*s3*h);
    e4=exp(k*s4*h);
    %A=[m1+n2 m2+n2 m3+n2 m4+n2 n2-n1; k1+k6 k2+k6 k3+k6 k4+k6 k5-k6; k7+k12 k8+k12 k9+k12 k10+k12 k11-k12; k1*e1 k2*e2 k3*e3 k4*e4 0; k7*e1 k8*e2 k9*e3 k10*e4 0];
    A1=[k1+k6 k2+k6 k3+k6 k4+k6; k7+k12 k8+k12 k9+k12 k10+k12; k1*e1 k2*e2 k3*e3 k4*e4; k7*e1 k8*e2 k9*e3 k10*e4];
    A2=[m1+n2 m2+n2 m3+n2 m4+n2; k7+k12 k8+k12 k9+k12 k10+k12; k1*e1 k2*e2 k3*e3 k4*e4; k7*e1 k8*e2 k9*e3 k10*e4];
    A3=[m1+n2 m2+n2 m3+n2 m4+n2; k1+k6 k2+k6 k3+k6 k4+k6; k1*e1 k2*e2 k3*e3 k4*e4; k7*e1 k8*e2 k9*e3 k10*e4];
    A1A1 =sym('A1A1',[4 4]);
    detA1=det(A1A1);
    D1=subs(detA1, A1A1, A1);
    A2A2 =sym('A2A2',[4 4]);
    detA2=det(A2A2);
    D2=subs(detA2, A2A2, A2);
    A3A3 =sym('A3A3',[4 4]);
    detA3=det(A3A3);
    D3=subs(detA3, A3A3, A3);
    D11=(n2-n1)*D1;
    D21=(k5-k6)*D2;
    D31=(k11-k12)*D3;
    D=D11-D21+D31;
    for k=4:0.2:5
        Dfun=matlabFunction(D);
        result=solve(Dfun==0,c);
        plot(k,real(result))
    end
end
hold on
for h=2.6
    e1=exp(k*s1*h);
    e2=exp(k*s2*h);
    e3=exp(k*s3*h);
    e4=exp(k*s4*h);
    %A=[m1+n2 m2+n2 m3+n2 m4+n2 n2-n1; k1+k6 k2+k6 k3+k6 k4+k6 k5-k6; k7+k12 k8+k12 k9+k12 k10+k12 k11-k12; k1*e1 k2*e2 k3*e3 k4*e4 0; k7*e1 k8*e2 k9*e3 k10*e4 0];
    A1=[k1+k6 k2+k6 k3+k6 k4+k6; k7+k12 k8+k12 k9+k12 k10+k12; k1*e1 k2*e2 k3*e3 k4*e4; k7*e1 k8*e2 k9*e3 k10*e4];
    A2=[m1+n2 m2+n2 m3+n2 m4+n2; k7+k12 k8+k12 k9+k12 k10+k12; k1*e1 k2*e2 k3*e3 k4*e4; k7*e1 k8*e2 k9*e3 k10*e4];
    A3=[m1+n2 m2+n2 m3+n2 m4+n2; k1+k6 k2+k6 k3+k6 k4+k6; k1*e1 k2*e2 k3*e3 k4*e4; k7*e1 k8*e2 k9*e3 k10*e4];
    A1A1 =sym('A1A1',[4 4]);
    detA1=det(A1A1);
    D1=subs(detA1, A1A1, A1);
    A2A2 =sym('A2A2',[4 4]);
    detA2=det(A2A2);
    D2=subs(detA2, A2A2, A2);
    A3A3 =sym('A3A3',[4 4]);
    detA3=det(A3A3);
    D3=subs(detA3, A3A3, A3);
    D11=(n2-n1)*D1;
    D21=(k5-k6)*D2;
    D31=(k11-k12)*D3;
    D=D11-D21+D31;
    for k=4:0.2:5
        Dfun=matlabFunction(D);
        result=solve(Dfun==0,c);
        plot(k,real(result))
    end
end
hold on
for h=2.8
    e1=exp(k*s1*h);
    e2=exp(k*s2*h);
    e3=exp(k*s3*h);
    e4=exp(k*s4*h);
    %A=[m1+n2 m2+n2 m3+n2 m4+n2 n2-n1; k1+k6 k2+k6 k3+k6 k4+k6 k5-k6; k7+k12 k8+k12 k9+k12 k10+k12 k11-k12; k1*e1 k2*e2 k3*e3 k4*e4 0; k7*e1 k8*e2 k9*e3 k10*e4 0];
    A1=[k1+k6 k2+k6 k3+k6 k4+k6; k7+k12 k8+k12 k9+k12 k10+k12; k1*e1 k2*e2 k3*e3 k4*e4; k7*e1 k8*e2 k9*e3 k10*e4];
    A2=[m1+n2 m2+n2 m3+n2 m4+n2; k7+k12 k8+k12 k9+k12 k10+k12; k1*e1 k2*e2 k3*e3 k4*e4; k7*e1 k8*e2 k9*e3 k10*e4];
    A3=[m1+n2 m2+n2 m3+n2 m4+n2; k1+k6 k2+k6 k3+k6 k4+k6; k1*e1 k2*e2 k3*e3 k4*e4; k7*e1 k8*e2 k9*e3 k10*e4];
    A1A1 =sym('A1A1',[4 4]);
    detA1=det(A1A1);
    D1=subs(detA1, A1A1, A1);
    A2A2 =sym('A2A2',[4 4]);
    detA2=det(A2A2);
    D2=subs(detA2, A2A2, A2);
    A3A3 =sym('A3A3',[4 4]);
    detA3=det(A3A3);
    D3=subs(detA3, A3A3, A3);
    D11=(n2-n1)*D1;
    D21=(k5-k6)*D2;
    D31=(k11-k12)*D3;
    D=D11-D21+D31;
    for k=4:0.2:5
        Dfun=matlabFunction(D);
        result=solve(Dfun==0,c);
        plot(k,real(result))
    end
end
hold off
toc

i try to draw the graph between real© and k for different value of h. but my programe is not run. please help

Answer 1

[Walter Roberson]

Replace

        Dfun=matlabFunction(D);
 

with

        Dfun = subs(D);

in multiple places.

[neetu malik]

Not work yet

[Walter Roberson]

Oh yes, that code. It is not possible to get that code to work.

I posted numeric code that does find solutions, several days ago, at https://www.mathworks.com/matlabcentral/answers/1773150-why-the-program-is-not-running#comment_2304010 .

There is NO chance that solve() will return a symbolic solution for you.

vpasolve() might potentially return a solution for you, if it were given a guess very close to the actual solution. But even given a very good starting point, it would take hours to find one solution.

Your current approach is not practical. The numeric approach, on the other hand, is practical.

Answer 2

[neetu malik]

This is also not give a proper solution of this problem

[Walter Roberson]

What is your numeric implementation of the current Question, rewritten along the lines of the numeric code that I posted in https://www.mathworks.com/matlabcentral/answers/1773150-why-the-program-is-not-running#comment_2304010 ?

Answer 3

[neetu malik]

sir i attached the numerical problem for which i want to make graph. please help to make graph of the problem.

[neetu malik]

i want the make figure 2

[Walter Roberson]

@neetu malik

Okay, go ahead and do that. But my advice from having worked with your equations is that you have no chance of doing this using solve() and no realistic chance of doing this using vpasolve() either. However if you write your code in numeric form the way I showed in https://www.mathworks.com/matlabcentral/answers/1773150-why-the-program-is-not-running#comment_2304010 then you have a chance.

The code I posted there already plots c against h and k in a 3D plot.