(* Set Constants *)
n = 1;
[Lambda] = 1;

(* Set functions *)
p[r_] = 1/r^2;
q[r_] = 1/r;

(* ODE's, which for some reason I called pde1/2 *)
pde1 = D[[Phi][r], {r, 2}] + q[r]*D[[Phi][r], {r, 1}] -
   p[r]*(n - [Alpha][r])^2 [Phi][r] + [Lambda]/
    2 (1 - [Phi][r]^2) [Phi][r] == 0
pde2 = D[[Alpha][r], {r, 2}] -
   q[r]*D[[Alpha][r], {r, 1}] + (n - [Alpha][r]) [Phi][r]^2 == 0

(* Min/max numbers *)
min = 0.000001;
max = 3.4; (* Define a "finite infinity" - I would like to set this higher *)

(* Boundary Conditions *)
bc = {[Phi][min] == 0, [Alpha][min] == 0, [Phi][max] ==
   1 - 1/max Exp[-max], [Alpha][max] == 1 - 5/max Exp[-max]};
(* solution[r_]=NDSolve[{pde1,pde2,bc},{[Phi],[Alpha]},{r,min,max},
Method[Rule]{StiffnessSwitching},AccuracyGoal[Rule]5,PrecisionGoal
[Rule]4] *)

sol = NDSolve[{pde1, pde2, bc}, {[Phi], [Alpha]}, {r, min, max},
  Method -> {"StiffnessSwitching",
    Method -> {"ExplicitRungeKutta", Automatic}}, AccuracyGoal -> 15,
  PrecisionGoal -> 15, MaxSteps -> 20000]

Plot[Evaluate[{[Phi][r], [Alpha][r]} /. sol], {r, 0, 4}]