See Carl's idea using  viewpoint  option  here

In this thread are also discussed some other ways to solve the problem.

There is substantial error in procedure  Simpson  - instead  n  should be  n/2 . I made a few less significant corrections. Now everything is working correctly:

Simpson := proc(f, a, b, n::even)

local h, S1, S2, S, i;

h := (b-a)/n;

S1 := 0.;

for i from 0 to n/2-1 do

S1 := S1 + f(a + (2*i+1)*h);

end do;

S2 := 0.;

for i from 1 to n/2-1 do

S2 := S2 + f(a + (2*i)*h);

end do;

evalf(h/3 * ( f(a)+f(b) + 4*S1 + 2*S2 )); 

end proc:

Example of work of the procedure  Simpson  and independent calculation:

Digits := 5;                          

f := x -> 1/sqrt(39.24*x-44.65*(x*arccos(x)-sqrt(1-x^2)-13.88*(1-x^2)^1.5));

Simpson(f, 0, 1, 100);

evalf(Int(f(x), x=0..1));

is(diff(abs(x), x) = x/abs(x))  assuming x < 0;

is(diff(abs(x), x) = x/abs(x))  assuming x > 0;

                                  true

                                  true

Slightly more programmatically:

op(-1, eval(splcurve) assuming v < 4.2);

If you want to get a normal manual solution, you can use  Student[Calculus1][ShowSolution]  command.

Example:

Student[Calculus1][ShowSolution](Int(sin(x)*(x + sin(x)),x));

roll:=rand(1..6):

L:=[seq([seq(roll(), j=1..10)], i=1..30)];

NumbersOfThrees:=map(t->ListTools[Occurrences](3, t), L);

The number  pi  should be coded as  Pi .

All curves are in the same plot:

restart;

Eq:=ln(2*s*t+2*s*x*sqrt(t)/sqrt(Pi)+x^2)-2*s*(arctan((s+x*sqrt(Pi/t))/sqrt(2*Pi*s-s^2))-(1/2)*Pi)/sqrt(2*Pi*s-s^2) = 0;

S:=[0.01,0.005,0.003,0.002];

plots[implicitplot]([seq(eval(Eq,s=S[i]), i=1..nops(S))], t=0..5, x=-2..2, color=[red,blue,green,yellow], gridrefine=5);

It is easy to verify identity

is(sin(2*x)^2 - cos(8*x)^2 = - cos(10*x)*cos(6*x));

                                         true

Therefore, the equation splits into two equations, which are easy to solve

solve({x >= 0, `or`(cos(10*x) = 0, cos(6*x)+1/2 = 0), x <= (1/2)*Pi}, AllSolutions, explicit);

Consider the function which is defined by the equation and its symmetries:

restart;

f:=x->sin(2*x)^2-cos(8*x)^2-(1/2)*cos(10*x):

is(f(Pi-x)=f(Pi+x)) and is(f(Pi/2-x)=f(Pi/2+x));

                                        true

Therefore it is sufficient to solve this equation in the range  0 .. Pi/2

RootFinding[Analytic](f(x), x, re = 0 .. (1/2)*Pi, im = -1 .. 1);

identify([%]);

You can enter a variable (the counter), first assign it the value  0 , and then at each step of the loop is increased by 1 .

Example: find all Pythagorean triangles with hypotenuse not exceeding 100 . Solution in the for loop:

restart;

k:=0: # The counter for the steps of the loop

n:=0: # The counter for solutions

for c from 2 to 100 do

for a from 1 to c-1 do

k:=k+1:

b:=sqrt(c^2-a^2):

if type(b,integer) and b>a then n:=n+1; L[n]:=[a,b,c] fi;

od:

od:

L:=convert(L, list):

k;  n;  L;

Addition:  The example of a purely illustrative nature. There are far more effective methods of generating Pythagorean triples.

restart;

 L:=[seq(seq(x^i*y^j, j=0..3-i), i=0..3)];

 add(a[k]*L[k+1], k=0..nops(L)-1);

The procedure  Euler  solves your problem:

restart;

Euler:=proc(Eq,Ivp,h,n)  # n is the number of steps

local Y, k;

Y[0]:=Ivp[2];

for k to n do

Y[k]:=Y[k-1]+h*eval(rhs(Eq),{x=Ivp[1]+h*(k-1), y(x)=Y[k-1]});

od;

[seq([Ivp[1]+i*h,Y[i]], i=0..n)];  # List of points of the approximate solution

end proc:

Example of use:

Eq:=diff(y(x),x)=(y(x)-x)^2; Ivp:=[0,0];

Euler(Eq, Ivp, 0.1, 10);

plot([%, x-tanh(x)], x=0..1, color=[red, blue]); # Graphs of approximate (red) and exact (blue) solutions

The code of the procedure was edited.

First get more awkward answer, but after a few conversions manages to simplify it.

restart;

dsolve({diff(y(x),x)=(y(x)-x)^2, y(0)=0}, y(x));

assign(%):

convert(expand(simplify(convert(y(x), trig))), tanh);

restart;

A:=(n,a,b,c,p)->((1+c*p)*a/(1-p)+(2+c*p)*b)*mul((1+c*(1-p)^k)*p, k=1..n-1)+add((a+b+k*b*(1-p))*mul((1+c*(1-p)^j)*p,j=k..n-1), k=2..n);

for n do

a[n]:=A(n,0.2,0.09,0.57,0.3);

if a[n]>=1  then Sol:=[[n-1,a[n-1]], [n,a[n]]]; break; fi;

od:

Sol;

You choose what  answer to you is more appropriate  n=6  or  n=7

The code is edited.

The simple procedure A  constructs the characteristic matrix with your properties of any size. For example displayed  A(10). Then we find the characteristic polynomial for  n = 50, and is set to 0 the coefficient of the imaginary part. We find that it is equal to  0  only for  x = 0. Thus we have a symmetrical real matrix that have only real roots:

A:=n->Matrix(n, {(1,1)=-I*x-lambda, (n,n)=-I*x-lambda, seq((i,i)=-lambda, i=2..n-1), seq((i,i-1)=-1, i=2..n), seq((i-1,i)=-1, i=2..n)}):

A(10);

B:=sort(LinearAlgebra[Determinant](A(50)), lambda);

x=solve(coeff(B,I)=0, x);

Verification:

C:=subs(x=0, B):

fsolve(C=0, lambda, complex);

We got only real roots.