@acer 

Hi again, I edited my equations and initial conditions, now i want to plot V(t) in terms of U0 in t=1s. tnx

eq1:=diff(r[1](t), t, t)+(.3754075596+.3045306134*U[0])*(diff(r[1](t), t))+740.0274933*r[1](t)-0.1181114312e-2*V(t) = -0.5234785518e-2*U[0]^2*q(t):
eq2:=2.030000000*10^(-8)*(diff(V(t), t))+4.065040650*10^(-7)*V(t)+0.1181114312e-2*(diff(r[1](t), t)) = 0:
eq3:=diff(q(t), t, t)-4.997988315*U[0]*(diff(q(t), t))+277.5543022*U[0]^2*q(t) = 0:
ics := V(0) = 0, q(0) = 0.1e-3, (D(q))(0) = 0, r[1](0) = 0, (D(r[1]))(0) = 0

@acer 

tnx for reply, forgive me for this wrong, the value A=1 and I just obtain r1(t), V(t) , q(t) in terms of U0, ignore the problem plot. I just need solve the ODEs in terms of U0. tnx. 

@mehdi jafari 

tnx for your answer, yes, your think is right, i have to obtain A_max where is the maximum value of the f(t) function at steady state. 

@Adri van der Meer 

thank you

@Christopher2222 

tnx

@gkokovidis 

tnx

@acer 

tnx my friend, i got it. 

@acer 

tnx my friend, i have one problem with ur methods, i try roots : 22.371016133870070862, 50.405610019621053775, 61.672103588567668205, but i can't earn lhs(eq)=0 with these roots.why?

for more explanation if we take x=22.371016133870070862 we earn lhs(eq)=8.2424325*10^8 whlie it must be zero. tnx again. 

@acer 

thank u so much... good luck.

@gaurav_rs 

tnx for ur answer my dear friend, but i must find the other roots. (assuming x>0)

@vv 

thank u so much my dear friend for ur excellent answer, can i solve too these equtions with ur method ?

eq1 := diff(u1(x), x, x)+2*x*(diff(u1(x), x))-u1(x)-(int(u1(s)-u2(s), s = 0 .. x)) = 2+x-exp(x)+2*x*exp(x)-cos(x);
eq2 := diff(u2(x), x, x)+diff(u2(x), x)-2*x*u2(x)-(int(u1(s)+u2(s), s = 0 .. x)) = 2*cos(x)-3*x-(1+2*x)*sin(x)-exp(x);
bcs := u1(0) = 1, (D(u1))(0) = 1, u2(0) = 1, (D(u2))(0) = 1

@Kitonum 

Thank u so much my friend. tnx.

@Preben Alsholm 

Hi, tnx for ur answer, the parameters M, R and aplha are constant, u can assume them M=R=alpha=1 or every value else, i know that must solve it with numerical method but how?

@Carl Love 

thanks my friend, You assume theta=0..Pi, but for theta=0..2Pi how? my plot for theta=0..2*Pi is like theta=0..Pi and this is not true. 

@mehran rajabi 

For more explanation my good friend, the range of Z for this problem is between -0.5 to 0.5 and I will write for loop so first it create values of Z between -0.5 to 0.5 and second the values of function f(z) in this interval([-0.5,0.5]) are closer together. I don't know what some of the values of Z create value closer together for function f(z) and I will obtain this values with maple.