Hi

However I run similar codes in Maple and Matlab, but the different results are observed.

By increasing number of basis (i.e. m) The convergence is observed, but results are different from each other.

The convergence for Maple: (m,P) =(10,10.8154),(20,10.8081),(30,10.8067),(50,10.8061)

The convergence for Matlab: (m,P) =(10,10.0474),(20,10.0194),(30,10.0143),(40,10.0125)

I have two questions:

1-which answer is correct?

2-why the different results for the same number of basis are obtained?

**Maple code (Maple 2016)**

restart;

tm := time():

with(LinearAlgebra):

Digits := 500:

beta := 1:

nu := 0.3:

lambda := 2:

G := 5:

ko := .5*0.1e7:

Ec := 0.380e12:

Em := 0.70e11:

ri := 0:

ro := 0.5:

ti := 0.1e-1:

K := 0.4e-1*0.1e10:

n := 1:

m := 20:

alpha := 0.1:

t := ti*(1+alpha*r/ro)^beta:

Er := (Ec-Em)*(r/ro)^n+Em:

W := simplify(add(a[n]*ChebyshevT(n, r), n = 0 .. m)):

sys := {eval(W, r = ro), eval(diff(W, r), r = ri)}:

W := subs(solve(sys, {a[0], a[1]}), W):

d1 := diff(W, r):

d2 := diff(d1, r):

W := subs(solve({eval(ko*d1-Ec*ti^3*(1+alpha)^(3*beta)*(d1*nu+d2*r)/(12*r*(-nu^2+1)), r = ro) = 0}, {a[2]}), W):

d1 := diff(W, r):

d2 := diff(d1, r):

Uf := int(K*(1+lambda*r/ro+G*(r/ro)^2)*W^2*r, r = ri .. ro):

NUM := int(((d2+d1/r)^2-(2*(1-nu))*d2*d1/r)*Er*t^3*r/(12*(-nu^2+1)), r = ri .. ro)+ko*ro*(eval(d1, r = ro))^2+Uf:

DEN := int(d1^2*r, r = ri .. ro):

PT := NUM-Em*ti^3*P*DEN/ro^2:

for c from 3 to m do

eq[c] := diff(PT, a[c]) end do:

MT := GenerateMatrix([seq(eq[j], j = 3 .. m)], [a[j]$j = 3 .. m])[1]:

solve(Determinant(MT)):

evalf[10](%[1]);

____________________________________________________________________________

**Matlab Code (R2014a)**

clear all;

digits(10);

tic

m = 20;

n = 1;

alpha = 0.1;

beta = 1;

nu = 0.3;

lambda = 2;

G = 5;

ko = 0.5*1e6;

Ec = 0.380e12;

Em = 0.70e11;

ro = 0.5;

ti = 0.1e-1;

K = 0.04*0.1e10;

ri=0;

syms r P;

c=sym('c%d',[1,m+1]);

c(1)=1;

c(2)=r;

for j=3:m+1

c(j)=simplify(2*r*c(j-1)-c(j-2));

end

Er = (Ec-Em)*(r/ro)^n+Em;

a=sym('a%d',[1,m+1]);

t = ti*(1+alpha*r/ro)^beta;

W = simplify(a*transpose(c));

[a(1),a(2)] = solve(subs(W,r,ro)==0,subs(diff(W,r),r,ri)==0,a(1),a(2));

W = simplify(a*transpose(c));

d1 = diff(W, r);

d2 = diff(d1, r);

a(3)=solve(subs(ko*d1-Ec*ti^3*(1+alpha)^(3*beta)*(d1*nu+d2*r)/(12*r*(-nu^2+1)),r,ro)==0,a(3));

W = simplify(a*transpose(c));

d1 = diff(W, r);

d2 = diff(d1, r);

Uf = int(K*(1+lambda*r/ro+G*(r/ro)^2)*W^2*r,r,ri,ro);

NUM = int(((d2+d1/r)^2-(2*(1-nu))*d2*d1/r)*Er*t^3*r/(12*(-nu^2+1)),r,ri,ro)+ko*ro*(subs(d1,r,ro))^2+Uf;

DEN = int(d1^2*r, r , ri , ro);

PT = NUM-Em*ti^3*P*DEN/ro^2;

for q=1:m-2

eq(q)= diff(PT, a(q+3));

end;

for q = 1:m-2

for u = 1:m-2

T = coeffs(eq(q),a(u+3));

tt(q,u) = T(2);

end;

end;

F=sort(vpa(solve(det(tt))));

F(1)

toc

___________________________________________________________