Hi Mark
Note that the homogeneous solution , phi_h, consists of harmonic functions, sines and cosines. Therefore their nth-derivitives are also harmonic but with an amplified magnitude with respect to thier arguments.
So, by substituting eq.7.100 in eq.7.99, with symbolic parameters and variables , Matlab executes 2nd and 4th derivitives of sine and cosine, which their magnitudes are in the form lambda^2 and lambda^4 respectively. In result , according to eq.100 , all terms will have a common factor of lambda^4, and so Matlab will represent just 2 terms of sine and cosine with lambda^4 factor and does not represent the form you want.
In fact, if you do the substitution manually, you'll get what you want, but Matlab applies symbolic derivation and shows the final result.
In your code, I suggest make all parameters symbolic, include, m, a, lambda_m, and then run your code.
I tried this, as below.
syms a m lamda_m x y phi_function_y_m phi_function_x_m
phi_h=phi_function_y_m*cos(lamda_m*x)+phi_function_x_m*sin(lamda_m*y);
deta_deta_phi_h=diff(diff(phi_h,x,2)+diff(phi_h,y,2),x,2)+diff(diff(phi_h,x,2)+diff(phi_h,y,2),y,2)
deta_deta_phi_h = 
You can also find the solution for phi_m numerically, according to the 4th-order ordinary differential equation for phi_m.
You could just write 2 mfiles, one of them is a function wherein you define your parameters and differential equation in the form of state-space, and the other one for simulation and solving the equation with a solver such 'ode45'.
Here is the state-space for the equation. After solving by ode solver, you could take the first element of state matrix, which is phi_m.