Basis for the Riemann-Roch space $L(kP)$ on a curve Let $C$ be a smooth algebraic plane complex curve of degree $d$ defined as the zero locus of a homogeneous polynomial $F$.
Let $P\in C$. For a positive integer $k$ consider the divisor 
$$D=kP$$ 

We know that the Riemann-Roch space $L(D)$ of meromorphic functions on $C$ with zeros prescribed and poles allowed by $D$ is a finite-dimensional complex vector space.

Question: Is there a standard way to find a basis of $L(kP)$ ?
 A: There's roughly two types of algorithms used to compute Riemann-Roch vector space bases - algebraic techniques and geometric techniques.  Assuming our curve is defined by an irreducible polynomial F(x,y) over a base field k, we wish to construct the vector space basis L(D) for some divisor D.
Algebraic Techniques
Assume the divisor has no points at infinity (change variables if it does).  Start with a basis for your entire algebraic function field as a $k(x)$-module, typically $1, y, \ldots, y^{n-1}$.  Multiply by powers of $(x-a)$, where $a$ lies under your divisor points, to put all of the basis functions into L(D).  Then compute series expansions at each of your divisor points, and use the coefficients from those expansions to form a certain matrix.  If its determinant is zero, then you can construct a linear combination that sums to zero and lets you reduce the order of a basis function while keeping it in L(D).  Keep going until the determinant is not zero.  Now you have a $k[x]$-basis for your divisor EXCEPT at infinity.  You look at your expansions at infinity to determine which powers of x you need to multiply by to get a k-basis for L(D).
Reference: Gilbert Ames Bliss, "Algebraic Functions", 1933, section 20.
KANT/KASH, MAGMA, and Sage all use an improved algorithm described in a 2001 paper by Florian Hess, "Computing Riemann-Roch spaces in algebraic function fields and related topics".  Hess's algorithm doesn't require the computation of series expansions, which can become a computational barrier.
Geometric Techniques
Geometric techniques are based on the observation that any $f \in L(D-E)$ induces a R-linear homomorphism from $L(-D)$ to $L(-E)$, where R is the coordinate ring of the curve, so $L(D-E) \cong Hom_R(L(-D), L(-E))$.  We can always pick D and E to be effective divisors, which allows us to identify $L(-D)$ and $L(-E)$ with ideals $I$ and $J$ in the coordinate ring, and now $L(D-E) \cong Hom_R(I, J)$.
This suggests that $L(D-E)$ should be closely related to the ideal quotient $J:I$.  We expect $L(D-E)$ to contain elements from the quotient field, however.  This can be remedied by picking a random element $f$ from $I$, computing $(fJ):I$, then dividing by $f$.  If the curve is non-singular, this basically works.  If the curve is singular, this approach needs to be modified by introducing an adjoint curve that intersects the original curve at its singularities.
Singular and Macaulay 2 both use the geometric approach involving adjoints.
Example
When I first answered this question in 2013, I did a sample calculation with KASH 2.5.  Now that Sage 9 has been released (Jan 2020), it makes sense to redo the example using Sage.
I'll compute basis elements for $L(nP)$ on the elliptic curve $y^2 = x^3-x$, where P is (0,0).
First, I set up the algebraic curve:
sage: R.<x> = FunctionField(QQbar)
sage: S.<Y> = R[]
sage: L.<y> = R.extension(Y^2 - (x^3 - x))

Even though the curve is specified using affine coordinates and an inhomogeneous polynomial, the calculation is done in projective space.
Now I verify that $(0,0)$ is an ordinary point of the curve:
sage: L.maximal_order().ideal(x,y).is_prime()
True

I use maximal_order to work with points in the finite portion of the plane.  There's also a maximal_order_infinite that would let me work with points at infinity, on an equal footing with finite points.  Like I said, the calculation is done in projective space, even though we use affine coordinates.
For a smooth curve, the ideal will always be prime.  In general, the ideal would not be prime if the point were a singularity.  In that case, I would use L's decomposition method to "resolve" the singularity, as they say, but in this example I can just keep going and form the divisor without having to resolve any singularities:
sage: D = L.maximal_order().ideal(x,y).divisor()

Now I can just use the Hess algorithm to compute basis elements for the Riemann-Roch spaces:
sage: D.basis_function_space()
[1]
sage: (2*D).basis_function_space()
[1, 1/x]
sage: (3*D).basis_function_space()
[1, 1/x, 1/x^2*y]
sage: (4*D).basis_function_space()
[1, 1/x, 1/x^2, 1/x^2*y]

The results are:
$$L(1P) = \{1\}\\
L(2P) = \{1, \frac{1}{x}\}\\
L(3P) = \{1, \frac{1}{x}, \frac{y}{x^2}\}\\
L(4P) = \{1, \frac{1}{x}, \frac{1}{x^2}, \frac{y}{x^2}\}$$
