'solve' seems to give me a solution that doesn't satisfy one of the equations in the system

I'm using the function 'solve' to solve a system of six equations and six variables. In the solution it gives, at least one of the equations seems not to be satisfied. Here's my code for the problem:

% given parameters
alpha = .05;        % exogenous rate of cost reductions
sigma = .05;         % exogenous uncertainty (infinitessimal st dev)
p0 = 10;             % initial carbon price
p1 = 50;            % carbon price under policy
lambda_p = .05;     % Poisson arrival rate of policy
lambda_q = .01;     % Probability of reversion to no policy state
r = .03;            % discount rate

% exponents for solution
beta1 = (.5*sigma^2 + alpha + sqrt((.5*sigma^2 + alpha)^2 + 2*sigma^2*(r+lambda_p)))/sigma^2;
beta2 = (.5*sigma^2 + alpha - sqrt((.5*sigma^2 + alpha)^2 + 2*sigma^2*(r+lambda_p)))/sigma^2;
beta3 = (.5*sigma^2 + alpha - sqrt((.5*sigma^2 + alpha)^2 + 2*sigma^2*r))/sigma^2;
beta4 = (.5*sigma^2 + alpha - sqrt((.5*sigma^2 + alpha)^2 + 2*sigma^2*(r+lambda_p+lambda_q)))/sigma^2;

syms a0 a1 ka ks B1 B2

S = solve(p1 - a1 == (lambda_p*lambda_q*ka*a1^beta3 + lambda_q*ks*a1^beta4)/(lambda_p + lambda_q), ...
    -1 == (lambda_p*lambda_q*ka*beta3*a1^(beta3-1) + lambda_q*ks*beta4*a1^(beta4-1))/(lambda_p + lambda_q), ... 
    p0 - a0 == B1*a0^beta1 + B2*a0^beta2 - lambda_p*a0/(alpha+lambda_p+r) + lambda_p*p1/(r+lambda_p), ... 
    -1 == B1*beta1*a0^(beta1-1) + B2*beta2*a0^(beta2-1) - lambda_p/(alpha+lambda_p+r), ... 
    (lambda_p*lambda_q*ka*a1^beta3 - lambda_p*ks*a1^beta4)/(lambda_p + lambda_q) == B1*a1^beta1 + B2*a1^beta2 - lambda_p*a1/(alpha+lambda_p+r) + lambda_p*p1/(r+lambda_p), ... 
    (lambda_p*lambda_q*ka*beta3*a1^(beta3-1) - lambda_p*ks*beta4*a1^(beta4-1))/(lambda_p + lambda_q) == B1*beta1*a1^(beta1-1) + B2*beta2*a1^(beta2-1) - lambda_p/(alpha+lambda_p+r));

The third equation isn't satisfied. When I obtain the output S, I calculate

>> p0-S.a0

ans =

-51.139888401602841212761952348352

and

>> S.B1*S.a0^beta1 + S.B2*S.a0^beta2 - lambda_p*S.a0/(alpha+lambda_p+r) + lambda_p*p1/(r+lambda_p)

ans =

4.58887467530944339251344448794

which aren't equal. Any ideas as to what's wrong? Thanks.