In the years since I left the answer below, Sage has improved dramatically, especially its symmetric function theory. In order to compute the third exterior power of the $S_7$-irrep corresponding to the partition $3+3+1$, one need only type the following:
s = SymmetricFunctions(QQ).schur()
s[1,1,1].inner_plethysm(s[3,3,1])
which returns
s[2, 1, 1, 1, 1, 1] + 2*s[2, 2, 1, 1, 1] + 4*s[2, 2, 2, 1] + 6*s[3, 1, 1, 1, 1] + 9*s[3, 2, 1, 1] + 7*s[3, 2, 2] + 4*s[3, 3, 1] + 7*s[4, 1, 1, 1] + 9*s[4, 2, 1] + 4*s[4, 3] + 2*s[5, 1, 1] + 4*s[5, 2] + s[6, 1] + s[7]
the same answer discovered by the GAP code below. Sage is also faster. From what I understand, GAP still wins when you want to study representations over a finite field. (But it's not really a contest since GAP is included in Sage!) From Sage, you may start a GAP console with
gap_console()
The answer from 2012 appears below.
Here is the GAP code I use to do these computations:
SchurFunctorOfCharacter:=function(char,p)
local n,t,c;
if p=[] then
return TrivialCharacter(UnderlyingCharacterTable(char));
fi;
n:=Sum(p);
t:=CharacterTable("Symmetric",n);
c:=List(CharacterParameters(t),u->u[2]);
return Symmetrizations([char],n)[Position(c,p)];
end;;
CharacterFromPartition:=function(table,p)
local c;
c:=List(CharacterParameters(table),u->u[2]);
return Irr(table)[Position(c,p)];
end;;
DecomposeCharacter:=char->List(Irr(UnderlyingCharacterTable(char)),x->ScalarProduct(x,char));;
t:=CharacterTable("Symmetric",7);;
chi:=CharacterFromPartition(t,[3,3,1]);;
DecomposeCharacter(SchurFunctorOfCharacter(chi,[1,1,1]));
This code computes the character table of $S_7$, finds the character corresponding to the partition $3+3+1$, and applies the Schur functor corresponding to the partition $1+1+1$ (otherwise known as $\wedge^3$). Here is the result:
[ 0, 1, 2, 4, 6, 9, 7, 4, 7, 9, 4, 2, 4, 1, 1 ]
These are the multiplicities of the irreducible constituents of our character. The ordering on partitions is lexicographic. For example, to determine the meaning of the $6$, just take the fifth partition of $7$:
Partitions(7)[5];
The output shows that the coefficient $6$ appears before the L-shaped partition $3+1+1+1+1$:
[3,1,1,1,1]