My result

Look at this switheart..

From this one,

Numerical solution of KdV equation

I need to combine with this one

Analytical higher order KdV bore

And here we go...this is my solution

The elevation of u versus x at t=25

And the coding is this longggggg

> restart;
with(plots):
n:=500;
dx:=0.1;dt:=0.001;
a:=1;b:=0;

for i from 1 to n/2 do
x[i]:=(i-(n/2))*dx;
f[i]:=a;
end do:
for i from (n/2)+1 to n do
x[i]:=(i-(n/2))*dx;
f[i]:=b;
end do:

g1:=listplot([seq([j*dt,f[j]],j=3..n-2)],thickness=2,color=red):

#FUNCTION 1
fn1:= proc()
global g, a, s, t, e, k, n:
a:=1:t:=1:e:=1:
e:=e+(g^2):
fn2():
while (evalf(t*((a^2)-s)) >= 6*(10^(-8))) do
fn2():
end do:
k:=evalf(Pi/(2*g)):
e:=(e*k)/2:
while (n <> 0) do
fn3():
end do:
n:=k:
g:=p:
end proc:

#FUNCTION 2
fn2:= proc()
global g, a, s, t, e, k, n:
s:=g*a:
a:=(a+g)/2:
g:=sqrt(s):
t:=2*t:
e:=evalf(e-t*((a^2)-s)):
end proc:

#FUNCTION 3
fn3:= proc()
global g, a, s, t, e, k, n:
n:=evalf(exp((-Pi*k)/n)):
end proc:

#MAIN FUNCTION
for i from 1 to 1000 do
n:=0:
p:=i/1000;
g:=evalf(sqrt(1-(p^2))):
fn1():
B:=1.0:
A:=0.0:
t1:=25:
C:=(B-A)*p+2*A-B+2*(B-A)*(e/k):
a:=(B-A)*p:
U[i]:=2*(B-A)*p+4*A+2*B-((4*(B-A)*p*(1-p)*k)/(e-(1-p)*k)):
K1[i]:=evalf((Pi/k)*((a/p)^(1/2))):
P[i]:=C+2*a*(((k-e)/p)/k):
R[i]:=C+2*a*((((k-e)/p)/k)-(1/p)):
Q[i]:=C+2*a*((((k-e)/p)/k)-1):
x[i]:=2*(B-A)*p+4*A+2*B-((4*(B-A)*p*(1-p))*k/(e-(1-p)*k)):
x1[i]:=x[i]/t1:
end do:

g2:=listplot([seq([x1[i],P[i]],i=1..999)],thickness=2,color=blue):
g3:=listplot([seq([x1[i],Q[i]],i=1..999)],thickness=2,color=green):

display(g1,g2,g3);

This is only 1% from all my works in this topic.. I wonder how long will it be if I need to describe the motion of the wave in future..=(