If you have Optimization Toolbox, you can use fsolve. Rewrite the three equations in the form F(x) = 0, where x is a vector of your three unknowns.
% set values
r0 = rand;
delta = rand;
sigma = rand;
P = rand;
% make function
f = @(x) [r0 + x(3)*delta + sigma*sqrt(delta) - x(1);
r0 + x(3)*delta - sigma*sqrt(delta) - x(2);
exp(-r0*delta)*(0.5*exp(-x(1)*delta) + 0.5*exp(-x(2)*delta))*100 - P];
% initial guess
x0 = rand(3,1);
% find solution
x0 = fsolve(f,x0)
Here, x(1) is r_u, x(2) is r_d, x(3) is theta.
EDIT TO ADD Having looked more closely at the math, I realized that you don't even need fsolve. The problem can be reduced to a single equation in one variable (theta):
% set values
r0 = rand;
delta = rand;
sigma = rand;
P = rand;
% make functions
r_u = @(theta) r0 + theta*delta + sigma*sqrt(delta);
r_d = @(theta) r0 + theta*delta - sigma*sqrt(delta);
f = @(theta) exp(-r0*delta)*(0.5*exp(-r_u(theta)*delta) + 0.5*exp(-r_u(theta)*delta))*100 - P;
% initial guess
theta0 = rand;
% find solution
theta0 = fzero(f,theta0)
r_u(theta0)
r_d(theta0)
