@Carl Love 

It didnt upload again. please see this code here, as you can see, it does not work. the last loop does not work ! 

restart;
> 
n:=20;
m:=1;
cc:=-200;
zzeta:=0.1;
sefr1:=0.3;
sefr:=0.2;
MM:=0;
lambda:=0.1;
Br:=1;
nn:=3;
> 
> 
> 
u(tau):=421.7129935*tau-2217.587728*tau^2+8897.376593*tau^3-27612.59182*tau^4+64248.00336*tau^5-1.083977605*10^5*tau^6-10.57029600-1.080951714*10^6*tau^13+7.999517316*10^5*tau^14-4.788741005*10^5*tau^15+2.309563748*10^5*tau^16+26511.11102*tau^18-5959.001794*tau^19+1.148523882*10^5*tau^7-95.23809524*tau^21+4.545454545*tau^22-9435.563781*tau^8-2.587683745*10^5*tau^9+6.473880128*10^5*tau^10+948.0272727*tau^20-88660.41892*tau^17-1.008692404*10^6*tau^11+1.175504242*10^6*tau^12;
> 
> 
B := 1+(2*(1-zzeta))*Br*(int(tau*(diff(u(tau), tau))^2, tau = (1/2)*zzeta*(1-zzeta) .. 1/2*(1-zzeta)));
eq4 := 4*B*u(tau)-(1+zzeta)*(diff(tau*(diff(theta(tau), tau)), tau))/tau-(1+zzeta)*Br*(diff(u(tau), tau))^2;


theta(tau):=sum(p^ii*theta[ii](tau),ii=0..nn);
HH:= p*((4*(1+(2*(1-zzeta))*Br*(int(tau*(diff(u(tau), tau))^2, tau = (1/2)*zzeta*(1-zzeta) .. 1/2-(1/2)*zzeta))))*u(tau)-(1+zzeta)*(diff(theta(tau), tau)+tau*(diff(theta(tau), tau, tau)))/tau-(1+zzeta)*Br*(diff(u(tau), tau))^2)+(1-p)*(diff(theta(tau),tau$2)):
eq5:=simplify(HH):
eq6:=collect(expand(eq5),p);

eq7:=
convert(series(collect(expand(eq5), p), p, nn+1), 'polynom');


for ii to nn do
ss[ii] := (coeff(eq7, p^ii)) ;
print (ii);
end do;

ss[0]:=diff(theta[0](tau), tau, tau);

icss[0]:=theta[0](0)=0, D(theta[0])(1)=1;

dsolve({ss[0], icss[0]});
theta[0](tau):= rhs(%);


for ii to nn do
ss[ii]:=(ss[ii]);
icss[ii]:=theta[ii](0)=0, D(theta[ii])(1)=0;
dsolve({ss[ii], icss[ii]});
theta[ii](tau):=value(rhs(%));
end do;





theta(tau):=evalf[10](sum((theta[jj])(tau),jj=0..nn));
> plot(, tau = .1 .. 1);

@Carl Love 

Dear carl
I attached the code, as you can see, it does not work 
please help me dama.mw

thanks again and again 

@Carl Love 

I dont know why I cant upload the picture, I copied and pasted the results here. 
thanks

for ii to nn do
theta[ii-1](tau):=convert(series(theta[ii-1](tau), tau, nn+1), 'polynom');
ss[ii]:=(ss[ii]);
icss[ii]:=theta[ii](0)=0, D(theta[ii])(1)=0;
dsolve({ss[ii], icss[ii]});
theta[ii](tau):=value(rhs(%));
end do;

tau
7 2 8 3
-5.297526110 10 tau + 3.892910316 10 tau

9 4 10 5
- 2.241308142 10 tau + 1.041425453 10 tau

10 6 13 13
- 3.969656447 10 tau - 1.144927761 10 tau

13 14 13 15
+ 2.558883005 10 tau - 4.182977586 10 tau

13 16 13 18
+ 4.803011179 10 tau - 5.366221145 10 tau

14 19 11 7
+ 1.948765755 10 tau + 1.245059033 10 tau

14 21 14 22
+ 5.916047990 10 tau - 7.541653420 10 tau

11 8 11 9
- 3.173675099 10 tau + 6.299381750 10 tau

11 10 14 20
- 8.396893740 10 tau - 3.872886396 10 tau

13 17 11 11
- 2.362225182 10 tau + 1.332932076 10 tau

12 12
+ 3.205721685 10 tau

/ d / d \\
+ 1.000000000 |----- |----- theta[1](tau)||
\ dtau \ dtau //

6 1.100000000 5
+ 5.364840735 10 tau - ----------- - 2.269580695 10
tau

14 23 14 24
+ 8.294399435 10 tau - 8.001553680 10 tau

14 25 14 26
+ 6.832053950 10 tau - 5.188740360 10 tau

14 27 14 29
+ 3.513543228 10 tau + 1.142903186 10 tau

13 30 12 32
- 5.473566470 10 tau - 8.699690700 10 tau

12 33 13 31
+ 2.854795441 10 tau + 2.323464346 10 tau

14 28 11 34
- 2.122529432 10 tau - 8.141338250 10 tau

11 35 10 36
+ 1.995748497 10 tau - 4.145120756 10 tau

9 37 8 38
+ 7.154081600 10 tau - 9.986090620 10 tau

8 39 6 40
+ 1.083350275 10 tau - 8.571320000 10 tau

5 41 42
+ 4.400000000 10 tau - 11000.00000 tau
theta[1](0) = 0, D(theta[1])(1) = 0
/ /tau
| | 453916139 1072968147 2
theta[1](tau) = | | --------- _z1 - ---------- _z1
| | 2000 400
\/0

529752611 3 973227579 4 2241308142 5
+ --------- _z1 - --------- _z1 + ---------- _z1
30 10 5

5207127265 6 3205721685000 13 5724638805000 14
- ---------- _z1 - ------------- _z1 + ------------- _z1
3 13 7

5117766010000 15 16
- ------------- _z1 + 2614360991250 _z1
3

3937041970000 18 53662211450000 19
+ ------------- _z1 + -------------- _z1
3 19

39696564470 7 129096213200000 21
+ ----------- _z1 + --------------- _z1
7 7

295802399500000 22 31126475825 8 317367509900 9
- --------------- _z1 - ----------- _z1 + ------------ _z1
11 2 9

10 20
- 62993817500 _z1 - 9743828775000 _z1

48030111790000 17 839689374000 11
- -------------- _z1 + ------------ _z1
17 11

12 754165342000000 23
- 11107767300 _z1 + --------------- _z1
23

103679992937500 24 25
- --------------- _z1 + 32006214720000 _z1
3

341602697500000 26 172958012000000 27
- --------------- _z1 + --------------- _z1
13 9

212252943200000 29 11429031860000 30
+ --------------- _z1 - -------------- _z1
29 3

32 2899896900000 33
- 726082608125 _z1 + ------------- _z1
11

54735664700000 31 87838580700000 28
+ -------------- _z1 - -------------- _z1
31 7

1427397720500 34 162826765000 35 36
- ------------- _z1 + ------------ _z1 - 5543745825 _z1
17 7

41451207560 37 3577040800 38 998609062 39
+ ----------- _z1 - ---------- _z1 + --------- _z1
37 19 39

43334011 40 8571320 41 220000 42 11000 43
- -------- _z1 + ------- _z1 - ------ _z1 + ----- _z1
16 41 21 43

\
11 | 340379736579462878110139
+ -- ln(_z1) d_z1| + ------------------------ tau
10 | 7849656799002018000
/

@Carl Love 

Dear carl
I think that maple shows the results because of the command " theta[ii-1](tau):=convert(series(theta[ii-1](tau), tau, nn+1), 'polynom'); " I want to get the results without it. As far as I know, this command forces the results to converge by eliminating some of them so the result is not accurate. I dont know why the image did not upload. I upload it again here.

thanks aloooooot for your time. :)

@Carl Love 

@Carl Love 

Thanks for your time. I really appreciate your favor. sadly it still does not work ! I attached the picture you understand better what I am saying. still some Z! appears.

@Carl Love 

Im waiting :)
thanks alot for your attention

@Carl Love 

The last dsolve in the code " dsolve({ss[ii], icss[ii]}); " connot integrate theta[1] twice and some "Z1" apears on the screen ! when I integrate theta [1] seprately, maple gives the best answer but it does not work in this code !

thanks for your attention

@Carl Love 


both will be helpfull, but the best opion would be Newton Raphson method

thanks

I applied this coomand but I did not get the best result. I know the values of "h1" and "bb" are in order of 0.01 but I have to obtain exact values of them. when I use fnormal just the second term of

will be omited from the equation. how can I omit every little values containing parameter of bb and h1 ?

I think if I can find a command that elemintaes every terms of the order of (10^-5) my problem will be solved

thanks

I applied this coomand but I did not get the best result. I know the values of "h1" and "bb" are in order of 0.01 but I have to obtain exact values of them. when I use fnormal just the second term of

will be omited from the equation. how can I omit every little values containing parameter of bb and h1 ?

I think if I can find a command that elemintaes every terms of the order of (10^-5) my problem will be solved

thanks

The difrential equation should be independent of all 3 constants, however, here there are 2 constants and jusc C is omitted. the result have to be thrd order diffrential equation.

Thanks for your attention

The difrential equation should be independent of all 3 constants, however, here there are 2 constants and jusc C is omitted. the result have to be thrd order diffrential equation.

Thanks for your attention