I was using pdsolve to solve for a 1-D longitudinal wave equation.  My particular problem added an external stimulus that is not mechanical in original, so the acoustic wave velocity has to be modeled as a variable c(H) and I also have to add a du/dx term in the model.  You can see my question posted on April 17, 2026.  I was told that pdsolve does not handle such a problem.  The suggestion in the posting is to use a finite difference method.  

I am verifying that approach by solving a simplier problem where c is constant but with an initial uniformly stretched material.  The solution does not seem to be physical.  The material should relax throughout the whole length of the material, but the solution shows relaxation at the end and stay uniformly stretched at the center.  I ran the problem with pdsolve and get a different result that I think is more realistic.

Is there something I can tweak in the finite difference approach to overcome that issue?

pde_finite_difference_method_linear_ic.mw

pdsolve_exercise_damping_ini_linear_a.mw

Please see the attached worksheet:

Question 1:  How to plot a partial portion of an array s2 using "dataplot," say the first 50 elements.

Question 2:  I have also tried to use "plot" to plot the first 50 elements s2 and keep on getting errors or strange looking plot with horizontal lines.  I have tried plot(s2,1..50), plot(s2,x=1..50), plot(s2,m=1..50), plot(s2[m],m=1..50), plot(s2[1..50],1..50], plot(1..50,s2[1..50]).

plot_exercise_a.mw

pdsolve is capable of numeric solution for a 1-D standing wave equation with a constant c:

diff(u(x, t), t, t) = c^2*diff(u(x, t), x, x)

In my problem, c is a function of lambda; both c and lambda has a common parameter H.  My problem has piecewise nonlinear function of c(H) and lambda(H) which cannot be reformulate as c(lambda) readily.  I therefore use a procedure based on fsolve to enter lambda and back out c.  In the example, I use piecewise linear function as demonstration.  I ran into a "(in fsolve) Can't handle expressions with typed procedures" error.

If I can overcome this, my input to the procedure also includes du/dx.  Can pdsolve handles such an equation?

I am not sure whether this should be in "create a post" or "ask a question." Let me know if this is more appropriate in "ask a question"

The standing wave equation is given by:

PDE:=diff(u(x,t),t,t)=c^2*diff(u(x,t),x,x)

IBC := {u(x, 0) = A0*cos(Pi*x/L), D[1](u)(0, t) = 0, D[1](u)(L, t) = 0, D[2](u)(x, 0) = 0}

pdsolve(PDE,IBC,numeric)

In my problem, the material is a magneto-elastic material where c, the speed of the acoustic wave, is a function of a magnetic field H.  The material is nonlinear and saturable.  I define it by a 3 segment piecewise nonlinear function of H.  The material response is a result of a sinusoidal H field.  I am interested in solving u(x,t).  

With that, c in the PDE has to be rewritten as c(H(t)), pdsolve gives an error as PDE has to be expressed as a function of u,t, or x.  So I redefine c(H(t)) as c(t).  I ran into another error in pdsolve as the piecewise has to be based on t or x and not H.  

The problem is that depending on H I can go through all 3 segments and back in one cycle and I have to find the corresponding t's for the piecewise.  Now I am driving the material with 100 cycles, I have to find and list all those piecewise transition points which is hardly practical.

Is there any other ways to approach and solve this problem?

These are the original set up to solve for second order partial differetial equations:

E := 1.456*10^11;
rho := 7900;
a := sqrt(E/rho);
L := 0.03639;

f0 := 1/(2*L)*a;
f := f0 - 50;
A0 := 11.75*L/2*10^(-6);

d := 0.8*1000;
T0 := 5/f;
NULL;
PDE := diff(u(x, t), t, t) = a^2*diff(u(x, t), x, x) - d*diff(u(x, t), t);
IBC := {u(x, 0) = A0*cos(Pi*x/L), D[1](u)(0, t) = 0, D[1](u)(L, t) = 0, D[2](u)(x, 0) = 0};
pds := pdsolve(PDE, IBC, numeric, timestep = 2*10^(-7));

 I know it solves properly because I can plot it using pds:-plot without problem.  The question is how to extract the value at specific x and t.  I started to try some of the more intuitive way such as u(0,0) without getting an answer.  I tried a number of eval expression such as eval(pds,[x=0,t=0]), eval(u(x,t),[x=0,t=0]) and so on without success.  

Then I tried using :-value as mentioned in an example from the internet, also resulting in an error.

val:=pds:-value(x=0,t=0)

It seems to be simple.  I guess I am too new with maple so be gentle.  Thank you.