I would keep it as the original system (and change ‘to’ to 0):
syms x(t) z(t) a0 a1 b w
Dz = diff(z);
D2z = diff(z,2);
Dx = diff(x);
D2x = diff(x,2);
Eq1 = D2x == -a0 * (a1 - b*Dz) * cos(w*t);
Eq2 = D2z == -a0 * b * Dx * cos(w*t);
Soln = dsolve(Eq1, Eq2, x(0) == 0, Dx(0) == 0, z(0) == 0, Dz(0) == 0);
X = Soln.x;
Z = Soln.z;
X = simplify(X, 'steps', 20)
Z = simplify(Z, 'steps', 20)
This gives you two ‘solutions’ involving integrals, that it transformed with dummy variable ‘y’:
X =
-(a1*int(sin((a0*b*sin(w*y))/w), y, 0, t))/b
Z =
-(a1*int(exp(-(a0*b*sin(w*y)*1i)/w)*(exp((a0*b*sin(w*y)*1i)/w) - 1)^2, y, 0, t))/(2*b)
This is likely as close as you can get to an analytic solution.
