(EDIT 1: Replaced hand-waving argument in third paragraph with a hopefully less incorrect version)

(EDIT 2: Added final paragraph about obtaining all conformal deformations for surfaces other than sphere.)

I think it is possible to see the infinitesimal rigidity of the sphere, even if it does involve a PDE as Dmitri says. I think you can also try and see if for other embedded surfaces, all infinitesimal deformations of conformal structure are accounted for by deformations of the embedding in a similar way.

For the case of S^{2}, what you want is to do is take a normal vector field V (i.e. infinitesimal change of embedding) and produce a tangent vector field X such that flowing along X gives the same infinitesimal change in conformal structure as flowing along V. This should amount to solving a linear PDE, so as Dmitri says a PDE is definitely involved, but probably not as hard as proving the existence of isothermal coordinates (which from memory is non-linear). For the standard embedding of S^{2} there can't be too many choices for this linear differential operator given that it has to respect the SO(3)-symmetry.

I guess we're looking for a first-order equivariant linear operator from normal vector fields to tangent vector fields. If we identify normal fields with functions then two possible candidates are to take X=grad V or X to be the Hamiltonian flow generated by V. I can't think of any others and probably it's possible to prove these are the only such ones. (Assuming it's elliptic, the symbol of the operator must be an SO(3)-equivariant isomorphism from T^{*}S^{2} to TS^{2} and there can't be too many choices! Using the metric leads to grad and using the area form leads to the Hamiltonian flow.) Then you just have to decide which one to use.

For the case of a general embedded surface $M$, you can ask "is it possible to obtain all deformations of conformal structure by deforming the embedding into R^{3}?" To answer this we can again think of a normal vector field as a function V on the surface. There is a second-order linear differential operator
$$
D\colon C^\infty(M) \to \Omega^{0,1}(T)
$$
which sends a normal vector field to the corresponding infinitesimal change of conformal structure. This operator will factor through the hessian with a homomorphism from $T^* \otimes T^*$ to $T^{*0,1}\otimes T^{1,0}$. The operator $D$ will not be onto, but what we want to know is whether every *cohomology class* in $H^{0,1}(T)$ has a representative in the image of $D$. At least, this is how I would try and approach the question; I'm sure there are other methods.