Here is a partial solution, but note that equation1_transformed and equation2_transformed are not satisfied for all values of t by the two solutions that the Symbolic Toolbox returns for the derived equations eqn1 and eqn2.
theta_val = 30*pi/180;
x_0_val = 10;
y_0_val = 20;
syms theta x_0 y_0 x_1 y_1 y__1 x__1 t x_2
equation1 = x_1 == x_0 + sin(theta + atan(x__1/y__1)) * sqrt(x__1^2 + y__1^2);
equation2 = y_1 == y_0 + cos(theta + atan(x__1/y__1)) * sqrt(x__1^2 + y__1^2);
equation1_sub = subs(equation1, [theta, x_0, y_0], [theta_val, x_0_val, y_0_val]);
equation2_sub = subs(equation2, [theta, x_0, y_0], [theta_val, x_0_val, y_0_val]);
ship_x = 2*t^2 + 2;
ship_y = 2*t + 3;
equation1_transformed = subs(equation1_sub, [x_1, y_1], [ship_x, ship_y])
equation2_transformed = subs(equation2_sub, [x_1, y_1], [ship_x, ship_y])
%solution = solve([equation1_transformed, equation2_transformed], [x__1, y__1]);
eqn1 = (2*t^2+2-10)^2+(2*t+3-20)^2 == x__1^2+y__1^2
eqn2 = tan(atan((2*t^2+2-10)/(2*t+3-20)) - pi/6) == x__1/y__1
solution = solve([eqn1,eqn2],[x__1,y__1])
solution.x__1
solution.y__1
simplify(subs([eqn1,eqn2],[x__1,y__1],[solution.x__1(1),solution.y__1(1)]))
simplify(subs([eqn1,eqn2],[x__1,y__1],[solution.x__1(2),solution.y__1(2)]))
simplify(subs([equation1_transformed,equation2_transformed],[x__1,y__1],[solution.x__1(1),solution.y__1(1)]))
simplify(subs([equation1_transformed,equation2_transformed],[x__1,y__1],[solution.x__1(2),solution.y__1(2)]))
